Context Navigation

mregular.lib @ 3686937

spielwiese

Last change on this file since 3686937 was 3686937, checked in by Oleksandr Motsak <motsak@…>, 11 years ago
Added '$Id$' as a comment to all libs (LIB/*.lib)
Property mode set to `100644`
File size: 57.0 KB

Line
1	//////////////////////////////////////////////////////////////////////////////
2	version="version mregular.lib 4.0.0.0 Jun_2013 "; // $Id$
3	category="Commutative Algebra";
4	info="
5	LIBRARY: mregular.lib Castelnuovo-Mumford regularity of homogeneous ideals
6	AUTHORS: I.Bermejo, ibermejo@ull.es
7	@* Ph.Gimenez, pgimenez@agt.uva.es
8	@* G.-M.Greuel, greuel@mathematik.uni-kl.de
9
10	OVERVIEW:
11	A library for computing the Castelnuovo-Mumford regularity of a homogeneous
12	ideal that DOES NOT require the computation of a minimal graded free
13	resolution of the ideal.
14	It also determines depth(basering/ideal) and satiety(ideal).
15	The procedures are based on 3 papers by Isabel Bermejo and Philippe Gimenez:
16	'On Castelnuovo-Mumford regularity of projective curves' Proc.Amer.Math.Soc.
17	128(5) (2000), 'Computing the Castelnuovo-Mumford regularity of some
18	subschemes of Pn using quotients of monomial ideals', Proceedings of
19	MEGA-2000, J. Pure Appl. Algebra 164 (2001), and 'Saturation and
20	Castelnuovo-Mumford regularity', Preprint (2004).
21
22	PROCEDURES:
23	regIdeal(id,[,e]); regularity of homogeneous ideal id
24	depthIdeal(id,[,e]); depth of S/id with S=basering, id homogeneous ideal
25	satiety(id,[,e]); saturation index of homogeneous ideal id
26	regMonCurve(li); regularity of projective monomial curve defined by li
27	NoetherPosition(id); Noether normalization of ideal id
28	is_NP(id); checks whether variables are in Noether position
29	is_nested(id); checks whether monomial ideal id is of nested type
30	";
31
32	LIB "general.lib";
33	LIB "algebra.lib";
34	LIB "sing.lib";
35	LIB "poly.lib";
36	//////////////////////////////////////////////////////////////////////////////
37	//
38	proc regIdeal (ideal i, list #)
39	"
40	USAGE: regIdeal (i[,e]); i ideal, e integer
41	RETURN: an integer, the Castelnuovo-Mumford regularity of i.
42	(returns -1 if i is not homogeneous)
43	ASSUME: i is a homogeneous ideal of the basering S=K[x(0)..x(n)].
44	e=0: (default)
45	If K is an infinite field, makes random changes of coordinates.
46	If K is a finite field, works over a transcendental extension.
47	e=1: Makes random changes of coordinates even when K is finite.
48	It works if it terminates, but may result in an infinite
49	loop. After 30 loops, a warning message is displayed and
50	-1 is returned.
51	NOTE: If printlevel > 0 (default = 0), additional info is displayed:
52	dim(S/i), depth(S/i) and end(H^(depth(S/i))(S/i)) are computed,
53	and an upper bound for the a-invariant of S/i is given.
54	The algorithm also determines whether the regularity is attained
55	or not at the last step of a minimal graded free resolution of i,
56	and if the answer is positive, the regularity of the Hilbert
57	function of S/i is given.
58	EXAMPLE: example regIdeal; shows some examples
59	"
60	{
61	//--------------------------- initialisation ---------------------------------
62	int e,ii,jj,H,h,d,time,lastv,sat,firstind;
63	int lastind,ch,nesttest,NPtest,nl,N,acc;
64	intmat ran;
65	def r0 = basering;
66	int n = nvars(r0)-1;
67	if ( size(#) > 0 )
68	{
69	e = #[1];
70	}
71	string s = "ring r1 = ",charstr(r0),",x(0..n),dp;";
72	execute(s);
73	ideal i,sbi,I,J,K,chcoord,m;
74	poly P;
75	map phi;
76	i = fetch(r0,i);
77	time=rtimer;
78	sbi=std(i);
79	ch=char(r1);
80	//----- Check ideal homogeneous
81	if ( homog(sbi) == 0 )
82	{
83	"// WARNING from proc regIdeal from lib mregular.lib:
84	// The ideal is not homogeneous!";
85	return (-1);
86	}
87	I=simplify(lead(sbi),1);
88	attrib(I,"isSB",1);
89	d=dim(I);
90	//----- If the ideal i is not proper:
91	if ( d == -1 )
92	{
93	dbprint(printlevel-voice+2,
94	"// The ideal i is (1)!
95	// Its Castelnuovo-Mumford regularity is:");
96	return (0);
97	}
98	//----- If the ideal i is 0:
99	if ( size(I) == 0 )
100	{
101	dbprint(printlevel-voice+2,
102	"// The ideal i is (0)!
103	// Its Castelnuovo-Mumford regularity is:");
104	return (0);
105	}
106	//----- When the ideal i is 0-dimensional:
107	if ( d == 0 )
108	{
109	H=maxdeg1(minbase(quotient(I,maxideal(1))))+1;
110	time=rtimer-time;
111	// Additional information:
112	dbprint(printlevel-voice+2,
113	"// Dimension of S/i : 0");
114	dbprint(printlevel-voice+2,
115	"// Time for computing regularity: " + string(time) + " sec.");
116	dbprint(printlevel-voice+2,
117	"// The Castelnuovo-Mumford regularity of i coincides with its satiety, and
118	// with the regularity of the Hilbert function of S/i. Its value is:");
119	return (H);
120	}
121	//----- Determine the situation: NT, or NP, or nothing.
122	//----- Choose the method depending on the situation, on the
123	//----- characteristic of the ground field, and on the option argument
124	//----- in order to get the mon. ideal of nested type associated to i
125	if ( e == 1 )
126	{ ch=0; }
127	NPtest=is_NP(I);
128	if ( NPtest == 1 )
129	{
130	nesttest=is_nested(I);
131	}
132	if ( ch != 0 )
133	{
134	if ( NPtest == 0 )
135	{
136	N=dn-d(d-1)/2;
137	s = "ring rtr = (ch,t(1..N)),x(0..n),dp;";
138	execute(s);
139	ideal chcoord,m,i,I;
140	poly P;
141	map phi;
142	i=imap(r1,i);
143	chcoord=select1(maxideal(1),1..(n-d+1));
144	acc=0;
145	for ( ii = 1; ii<=d; ii++ )
146	{
147	matrix trex[1][n-d+ii+1]=t((1+acc)..(n-d+ii+acc)),1;
148	m=select1(maxideal(1),1..(n-d+1+ii));
149	for ( jj = 1; jj<=n-d+ii+1; jj++ )
150	{
151	P=P+trex[1,jj]*m[jj];
152	}
153	chcoord[n-d+1+ii]=P;
154	P=0;
155	acc=acc+n-d+ii;
156	kill trex;
157	}
158	phi=rtr,chcoord;
159	I=simplify(lead(std(phi(i))),1);
160	setring r1;
161	I=imap(rtr,I);
162	attrib(I,"isSB",1);
163	}
164	else
165	{
166	if ( nesttest == 0 )
167	{
168	N=d*(d-1)/2;
169	s = "ring rtr = (ch,t(1..N)),x(0..n),dp;";
170	execute(s);
171	ideal chcoord,m,i,I;
172	poly P;
173	map phi;
174	i=imap(r1,i);
175	chcoord=select1(maxideal(1),1..(n-d+2));
176	acc=0;
177	for ( ii = 1; ii<=d-1; ii++ )
178	{
179	matrix trex[1][ii+1]=t((1+acc)..(ii+acc)),1;
180	m=select1(maxideal(1),(n-d+2)..(n-d+2+ii));
181	for ( jj = 1; jj<=ii+1; jj++ )
182	{
183	P=P+trex[1,jj]*m[jj];
184	}
185	chcoord[n-d+2+ii]=P;
186	P=0;
187	acc=acc+ii;
188	kill trex;
189	}
190	phi=rtr,chcoord;
191	I=simplify(lead(std(phi(i))),1);
192	setring r1;
193	I=imap(rtr,I);
194	attrib(I,"isSB",1);
195	}
196	}
197	}
198	else
199	{
200	if ( NPtest == 0 )
201	{
202	while ( nl < 30 )
203	{
204	chcoord=select1(maxideal(1),1..(n-d+1));
205	nl=nl+1;
206	for ( ii = 1; ii<=d; ii++ )
207	{
208	ran=random(100,1,n-d+ii);
209	ran=intmat(ran,1,n-d+ii+1);
210	ran[1,n-d+ii+1]=1;
211	m=select1(maxideal(1),1..(n-d+1+ii));
212	for ( jj = 1; jj<=n-d+ii+1; jj++ )
213	{
214	P=P+ran[1,jj]*m[jj];
215	}
216	chcoord[n-d+1+ii]=P;
217	P=0;
218	}
219	phi=r1,chcoord;
220	dbprint(printlevel-voice+2,"// (1 random change of coord.)");
221	I=simplify(lead(std(phi(i))),1);
222	attrib(I,"isSB",1);
223	NPtest=is_NP(I);
224	if ( NPtest == 1 )
225	{
226	break;
227	}
228	}
229	if ( NPtest == 0 )
230	{
231	"// WARNING from proc regIdeal from lib mregular.lib:
232	// The procedure has entered in 30 loops and could not put the variables
233	// in Noether position: in your example the method using random changes
234	// of coordinates may enter an infinite loop when the field is finite.
235	// Try removing this optional argument.";
236	return (-1);
237	}
238	i=phi(i);
239	nesttest=is_nested(I);
240	}
241	if ( nesttest == 0 )
242	{
243	while ( nl < 30 )
244	{
245	chcoord=select1(maxideal(1),1..(n-d+2));
246	nl=nl+1;
247	for ( ii = 1; ii<=d-1; ii++ )
248	{
249	ran=random(100,1,ii);
250	ran=intmat(ran,1,ii+1);
251	ran[1,ii+1]=1;
252	m=select1(maxideal(1),(n-d+2)..(n-d+2+ii));
253	for ( jj = 1; jj<=ii+1; jj++ )
254	{
255	P=P+ran[1,jj]*m[jj];
256	}
257	chcoord[n-d+2+ii]=P;
258	P=0;
259	}
260	phi=r1,chcoord;
261	dbprint(printlevel-voice+2,"// (1 random change of coord.)");
262	I=simplify(lead(std(phi(i))),1);
263	attrib(I,"isSB",1);
264	nesttest=is_nested(I);
265	if ( nesttest == 1 )
266	{
267	break;
268	}
269	}
270	if ( nesttest == 0 )
271	{
272	"// WARNING from proc regIdeal from lib mregular.lib:
273	// The procedure has entered in 30 loops and could not find a monomial
274	// ideal of nested type with the same regularity as your ideal: in your
275	// example the method using random changes of coordinates may enter an
276	// infinite loop when the field is finite.
277	// Try removing this optional argument.";
278	return (-1);
279	}
280	}
281	}
282	//
283	// At this stage, we have obtained a monomial ideal I of nested type
284	// such that reg(i)=reg(I). We now compute reg(I).
285	//
286	//----- When S/i is Cohen-Macaulay:
287	for ( ii = n-d+2; ii <= n+1; ii++ )
288	{
289	K=K+select(I,ii);
290	}
291	if ( size(K) == 0 )
292	{
293	s="ring nr = ",charstr(r0),",x(0..n-d),dp;";
294	execute(s);
295	ideal I;
296	I = imap(r1,I);
297	H=maxdeg1(minbase(quotient(I,maxideal(1))))+1;
298	time=rtimer-time;
299	// Additional information:
300	dbprint(printlevel-voice+2,
301	"// S/i is Cohen-Macaulay");
302	dbprint(printlevel-voice+2,
303	"// Dimension of S/i ( = depth(S/i) ): "+string(d));
304	dbprint(printlevel-voice+2,
305	"// Regularity attained at the last step of m.g.f.r. of i: YES");
306	dbprint(printlevel-voice+2,
307	"// Regularity of the Hilbert function of S/i: " + string(H-d));
308	dbprint(printlevel-voice+2,
309	"// Time for computing regularity: " + string(time) + " sec.");
310	dbprint(printlevel-voice+2,
311	"// The Castelnuovo-Mumford regularity of i is:");
312	return(H);
313	}
314	//----- When d=1:
315	if ( d == 1 )
316	{
317	H=maxdeg1(simplify(reduce(quotient(I,maxideal(1)),I),2))+1;
318	sat=H;
319	J=subst(I,x(n),1);
320	s = "ring nr = ",charstr(r0),",x(0..n-1),dp;";
321	execute(s);
322	ideal J=imap(r1,J);
323	attrib(J,"isSB",1);
324	h=maxdeg1(minbase(quotient(J,maxideal(1))))+1;
325	time=rtimer-time;
326	if ( h > H )
327	{
328	H=h;
329	}
330	// Additional information:
331	dbprint(printlevel-voice+2,
332	"// Dimension of S/i: 1");
333	dbprint(printlevel-voice+2,
334	"// Depth of S/i: 0");
335	dbprint(printlevel-voice+2,
336	"// Satiety of i: "+string(sat));
337	dbprint(printlevel-voice+2,
338	"// Upper bound for the a-invariant of S/i: end(H^1(S/i)) <= "+
339	string(h-2));
340	if ( H == sat )
341	{
342	dbprint(printlevel-voice+2,
343	"// Regularity attained at the last step of m.g.f.r. of i: YES");
344	dbprint(printlevel-voice+2,
345	"// Regularity of the Hilbert function of S/i: "+string(H));
346	}
347	else
348	{
349	dbprint(printlevel-voice+2,
350	"// Regularity attained at the last step of m.g.f.r. of i: NO");
351	}
352	dbprint(printlevel-voice+2,
353	"// Time for computing regularity: "+ string(time) + " sec.");
354	dbprint(printlevel-voice+2,
355	"// The Castelnuovo-Mumford regularity of i is:");
356	return(H);
357	}
358	//----- Now d>1 and S/i is not Cohen-Macaulay:
359	//
360	//----- First, determine the last variable really occuring
361	lastv=n-d;
362	h=n;
363	while ( lastv == n-d and h > n-d )
364	{
365	K=select(I,h+1);
366	if ( size(K) == 0 )
367	{
368	h=h-1;
369	}
370	else
371	{
372	lastv=h;
373	}
374	}
375	//----- and compute Castelnuovo-Mumford regularity:
376	s = "ring nr = ",charstr(r0),",x(0..lastv),dp;";
377	execute(s);
378	ideal I,K,KK,LL;
379	I=imap(r1,I);
380	attrib(I,"isSB",1);
381	K=simplify(reduce(quotient(I,maxideal(1)),I),2);
382	H=maxdeg1(K)+1;
383	firstind=H;
384	KK=minbase(subst(I,x(lastv),1));
385	for ( ii = n-lastv; ii<=d-2; ii++ )
386	{
387	LL=minbase(subst(I,x(n-ii-1),1));
388	attrib(LL,"isSB",1);
389	s = "ring mr = ",charstr(r0),",x(0..n-ii-1),dp;";
390	execute(s);
391	ideal K,KK;
392	KK=imap(nr,KK);
393	attrib(KK,"isSB",1);
394	K=simplify(reduce(quotient(KK,maxideal(1)),KK),2);
395	h=maxdeg1(K)+1;
396	if ( h > H )
397	{
398	H=h;
399	}
400	setring nr;
401	kill mr;
402	KK=LL;
403	}
404	// We must determine one more sat. index:
405	s = "ring mr = ",charstr(r0),",x(0..n-d),dp;";
406	execute(s);
407	ideal KK,K;
408	KK=imap(nr,KK);
409	attrib(KK,"isSB",1);
410	K=simplify(reduce(quotient(KK,maxideal(1)),KK),2);
411	h=maxdeg1(K)+1;
412	lastind=h;
413	if ( h > H )
414	{
415	H=h;
416	}
417	setring nr;
418	kill mr;
419	time=rtimer-time;
420	// Additional information:
421	dbprint(printlevel-voice+2,
422	"// Dimension of S/i: "+string(d));
423	dbprint(printlevel-voice+2,
424	"// Depth of S/i: "+string(n-lastv));
425	dbprint(printlevel-voice+2,
426	"// end(H^"+string(n-lastv)+"(S/i)) = "
427	+string(firstind-n+lastv-1));
428	dbprint(printlevel-voice+2,
429	"// Upper bound for the a-invariant of S/i: end(H^"
430	+string(d)+"(S/i)) <= "+string(lastind-d-1));
431	if ( H == firstind )
432	{
433	dbprint(printlevel-voice+2,
434	"// Regularity attained at the last step of m.g.f.r. of i: YES");
435	dbprint(printlevel-voice+2,
436	"// Regularity of the Hilbert function of S/i: "
437	+string(H-n+lastv));
438	}
439	else
440	{
441	dbprint(printlevel-voice+2,
442	"// Regularity attained at the last step of m.g.f.r. of i: NO");
443	}
444	dbprint(printlevel-voice+2,
445	"// Time for computing regularity: "+ string(time) + " sec.");
446	dbprint(printlevel-voice+2,
447	"// The Castelnuovo-Mumford regularity of i is:");
448	return(H);
449	}
450	example
451	{ "EXAMPLE:"; echo = 2;
452	ring r=0,(x,y,z,t,w),dp;
453	ideal i=y2t,x2y-x2z+yt2,x2y2,xyztw,x3z2,y5+xz3w-x2zw2,x7-yt2w4;
454	regIdeal(i);
455	regIdeal(lead(std(i)));
456	// Additional information is displayed if you change printlevel (=1);
457	}
458	////////////////////////////////////////////////////////////////////////////////
459	/*
460	Out-commented examples:
461	//
462	ring s=0,x(0..5),dp;
463	ideal i=x(2)^2-x(4)x(5),x(1)x(2)-x(0)x(5),x(0)x(2)-x(1)*x(4),
464	x(1)^2-x(3)x(5),x(0)x(1)-x(2)x(3),x(0)^2-x(3)x(4);
465	regIdeal(i);
466	// Our procedure works when a min. graded free resol. can
467	// not be computed. In this easy example, regularity can also
468	// be obtained using a m.g.f.r.:
469	nrows(betti(mres(i,0)));
470	ring r1=0,(x,y,z,t),dp;
471	// Ex.2.5 in [Bermejo-Gimenez], Proc.Amer.Math.Soc. 128(5):
472	ideal i = x17y14-y31, x20y13, x60-y36z24-x20z20t20;
473	regIdeal(i);
474	// Ex.2.9 in [Bermejo-Gimenez], Proc.Amer.Math.Soc. 128(5):
475	int k=43;
476	ideal j=x17y14-y31,x20y13,x60-y36z24-x20z20t20,y41z^k-y40z^(k+1);
477	regIdeal(j);
478	k=14;
479	j=x17y14-y31,x20y13,x60-y36z24-x20z20t20,y41z^k-y40z^(k+1);
480	regIdeal(j);
481	k=22;
482	j=x17y14-y31,x20y13,x60-y36z24-x20z20t20,y41z^k-y40z^(k+1);
483	regIdeal(j);
484	k=315;
485	j=x17y14-y31,x20y13,x60-y36z24-x20z20t20,y41z^k-y40z^(k+1);
486	regIdeal(j);
487	// Example in Rk.2.10 in [Bermejo-Gimenez], ProcAMS 128(5):
488	ideal h=x2-3xy+5xt,xy-3y2+5yt,xz-3yz,2xt-yt,y2-yz-2yt;
489	regIdeal(h);
490	// The initial ideal is not saturated
491	regIdeal(lead(std(h)));
492	// More examples:
493	i=y4-t3z, x3t-y2z2, x3y2-t2z3, x6-tz5;
494	regIdeal(i);
495	//
496	regIdeal(maxideal(4));
497	//
498	ring r2=0,(x,y,z,t,w),dp;
499	ideal i = xy-zw,x3-yw2,x2z-y2w,y3-xz2,-y2z3+xw4+tw4+w5,-yz4+x2w3+xtw3+xw4,
500	-z5+x2tw2+x2w3+yw4;
501	regIdeal(i);
502	//
503	ring r3=0,(x,y,z,t,w,u),dp;
504	ideal i=imap(r2,i);
505	regIdeal(i);
506	// Next example is the defining ideal of the 2nd. Veronesean of P3, a variety
507	// in P8 which is arithmetically Cohen-Macaulay:
508	ring r4=0,(a,b,c,d,x(0..9)),dp;
509	ideal i= x(0)-ab,x(1)-ac,x(2)-ad,x(3)-bc,x(4)-bd,x(5)-cd,
510	x(6)-a2,x(7)-b2,x(8)-c2,x(9)-d2;
511	ideal ei=eliminate(i,abcd);
512	ring r5=0,x(0..9),dp;
513	ideal i=imap(r4,ei);
514	regIdeal(i);
515	// Here is an example where the computation of a m.g.f.r. of I costs:
516	ring r8=0,(x,y,z,t,u,a,b),dp;
517	ideal i=u-b40,t-a40,x-a23b17,y-a22b18+ab39,z-a25b15;
518	ideal ei=eliminate(i,ab); // It takes a few seconds to compute the ideal
519	ring r9=0,(x,y,z,t,u),dp;
520	ideal i=imap(r8,ei);
521	regIdeal(i); // This is very fast.
522	// Now you can use mres(i,0) to compute a m.g.f.r. of the ideal!
523	//
524	// The computation of the m.g.f.r. of the following example did not succeed
525	// using the command mres:
526	ring r10=0,(x(0..8),s,t),dp;
527	ideal i=x(0)-st24,x(1)-s2t23,x(2)-s3t22,x(3)-s9t16,x(4)-s11t14,x(5)-s18t7,
528	x(6)-s24t,x(7)-t25,x(8)-s25;
529	ideal ei=eliminate(i,st);
530	ring r11=0,x(0..8),dp;
531	ideal i=imap(r10,ei);
532	regIdeal(i);
533	// More examples where not even sres works:
534	// Be careful: elimination takes some time here, but it succeeds!
535	ring r12=0,(s,t,u,x(0..14)),dp;
536	ideal i=x(0)-st6u8,x(1)-s5t3u7,x(2)-t11u4,x(3)-s9t4u2,x(4)-s2t7u6,x(5)-s7t7u,
537	x(6)-s10t5,x(7)-s4t6u5,x(8)-s13tu,x(9)-s14u,x(10)-st2u12,x(11)-s3t9u3,
538	x(12)-s15,x(13)-t15,x(14)-u15;
539	ideal ei=eliminate(i,stu);
540	size(ei);
541	ring r13=0,x(0..14),dp;
542	ideal i=imap(r12,ei);
543	size(i);
544	regIdeal(i);
545	*/
546	///////////////////////////////////////////////////////////////////////////////
547	///////////////////////////////////////////////////////////////////////////////
548	///////////////////////////////////////////////////////////////////////////////
549
550	proc depthIdeal (ideal i, list #)
551	"
552	USAGE: depthIdeal (i[,e]); i ideal, e integer
553	RETURN: an integer, the depth of S/i where S=K[x(0)..x(n)] is the basering.
554	(returns -1 if i is not homogeneous or if i=(1))
555	ASSUME: i is a proper homogeneous ideal.
556	e=0: (default)
557	If K is an infinite field, makes random changes of coordinates.
558	If K is a finite field, works over a transcendental extension.
559	e=1: Makes random changes of coordinates even when K is finite.
560	It works if it terminates, but may result in an infinite
561	loop. After 30 loops, a warning message is displayed and
562	-1 is returned.
563	NOTE: If printlevel > 0 (default = 0), dim(S/i) is also displayed.
564	EXAMPLE: example depthIdeal; shows some examples
565	"
566	{
567	//--------------------------- initialisation ---------------------------------
568	int e,ii,jj,h,d,time,lastv,ch,nesttest,NPtest,nl,N,acc;
569	intmat ran;
570	def r0 = basering;
571	int n = nvars(r0)-1;
572	if ( size(#) > 0 )
573	{
574	e = #[1];
575	}
576	string s = "ring r1 = ",charstr(r0),",x(0..n),dp;";
577	execute(s);
578	ideal i,sbi,I,J,K,chcoord,m;
579	poly P;
580	map phi;
581	i = fetch(r0,i);
582	time=rtimer;
583	sbi=std(i);
584	ch=char(r1);
585	//----- Check ideal homogeneous
586	if ( homog(sbi) == 0 )
587	{
588	"// WARNING from proc depthIdeal from lib mregular.lib:
589	// The ideal is not homogeneous!";
590	return (-1);
591	}
592	I=simplify(lead(sbi),1);
593	attrib(I,"isSB",1);
594	d=dim(I);
595	//----- If the ideal i is not proper:
596	if ( d == -1 )
597	{
598	"// WARNING from proc depthIdeal from lib mregular.lib:
599	// The ideal i is (1)!";
600	return (-1);
601	}
602	//----- If the ideal i is 0:
603	if ( size(I) == 0 )
604	{
605	dbprint(printlevel-voice+2,
606	"// The ideal i is (0)!
607	// The depth of S/i is:");
608	return (d);
609	}
610	//----- When the ideal i is 0-dimensional:
611	if ( d == 0 )
612	{
613	time=rtimer-time;
614	// Additional information:
615	dbprint(printlevel-voice+2,
616	"// Dimension of S/i : 0 (S/i is Cohen-Macaulay)");
617	dbprint(printlevel-voice+2,
618	"// Time for computing the depth: " + string(time) + " sec.");
619	dbprint(printlevel-voice+2,
620	"// The depth of S/i is:");
621	return (0);
622	}
623	//----- Determine the situation: NT, or NP, or nothing.
624	//----- Choose the method depending on the situation, on the
625	//----- characteristic of the ground field, and on the option argument
626	//----- in order to get the mon. ideal of nested type associated to i
627	if ( e == 1 )
628	{
629	ch=0;
630	}
631	NPtest=is_NP(I);
632	if ( NPtest == 1 )
633	{
634	nesttest=is_nested(I);
635	}
636	if ( ch != 0 )
637	{
638	if ( NPtest == 0 )
639	{
640	N=dn-d(d-1)/2;
641	s = "ring rtr = (ch,t(1..N)),x(0..n),dp;";
642	execute(s);
643	ideal chcoord,m,i,I;
644	poly P;
645	map phi;
646	i=imap(r1,i);
647	chcoord=select1(maxideal(1),1..(n-d+1));
648	acc=0;
649	for ( ii = 1; ii<=d; ii++ )
650	{
651	matrix trex[1][n-d+ii+1]=t((1+acc)..(n-d+ii+acc)),1;
652	m=select1(maxideal(1),1..(n-d+1+ii));
653	for ( jj = 1; jj<=n-d+ii+1; jj++ )
654	{
655	P=P+trex[1,jj]*m[jj];
656	}
657	chcoord[n-d+1+ii]=P;
658	P=0;
659	acc=acc+n-d+ii;
660	kill trex;
661	}
662	phi=rtr,chcoord;
663	I=simplify(lead(std(phi(i))),1);
664	setring r1;
665	I=imap(rtr,I);
666	attrib(I,"isSB",1);
667	}
668	else
669	{
670	if ( nesttest == 0 )
671	{
672	N=d*(d-1)/2;
673	s = "ring rtr = (ch,t(1..N)),x(0..n),dp;";
674	execute(s);
675	ideal chcoord,m,i,I;
676	poly P;
677	map phi;
678	i=imap(r1,i);
679	chcoord=select1(maxideal(1),1..(n-d+2));
680	acc=0;
681	for ( ii = 1; ii<=d-1; ii++ )
682	{
683	matrix trex[1][ii+1]=t((1+acc)..(ii+acc)),1;
684	m=select1(maxideal(1),(n-d+2)..(n-d+2+ii));
685	for ( jj = 1; jj<=ii+1; jj++ )
686	{
687	P=P+trex[1,jj]*m[jj];
688	}
689	chcoord[n-d+2+ii]=P;
690	P=0;
691	acc=acc+ii;
692	kill trex;
693	}
694	phi=rtr,chcoord;
695	I=simplify(lead(std(phi(i))),1);
696	setring r1;
697	I=imap(rtr,I);
698	attrib(I,"isSB",1);
699	}
700	}
701	}
702	else
703	{
704	if ( NPtest == 0 )
705	{
706	while ( nl < 30 )
707	{
708	chcoord=select1(maxideal(1),1..(n-d+1));
709	nl=nl+1;
710	for ( ii = 1; ii<=d; ii++ )
711	{
712	ran=random(100,1,n-d+ii);
713	ran=intmat(ran,1,n-d+ii+1);
714	ran[1,n-d+ii+1]=1;
715	m=select1(maxideal(1),1..(n-d+1+ii));
716	for ( jj = 1; jj<=n-d+ii+1; jj++ )
717	{
718	P=P+ran[1,jj]*m[jj];
719	}
720	chcoord[n-d+1+ii]=P;
721	P=0;
722	}
723	phi=r1,chcoord;
724	dbprint(printlevel-voice+2,"// (1 random change of coord.)");
725	I=simplify(lead(std(phi(i))),1);
726	attrib(I,"isSB",1);
727	NPtest=is_NP(I);
728	if ( NPtest == 1 )
729	{
730	break;
731	}
732	}
733	if ( NPtest == 0 )
734	{
735	"// WARNING from proc depthIdeal from lib mregular.lib:
736	// The procedure has entered in 30 loops and could not put the variables
737	// in Noether position: in your example the method using random changes
738	// of coordinates may enter an infinite loop when the field is finite.
739	// Try removing this optional argument.";
740	return (-1);
741	}
742	i=phi(i);
743	nesttest=is_nested(I);
744	}
745	if ( nesttest == 0 )
746	{
747	while ( nl < 30 )
748	{
749	chcoord=select1(maxideal(1),1..(n-d+2));
750	nl=nl+1;
751	for ( ii = 1; ii<=d-1; ii++ )
752	{
753	ran=random(100,1,ii);
754	ran=intmat(ran,1,ii+1);
755	ran[1,ii+1]=1;
756	m=select1(maxideal(1),(n-d+2)..(n-d+2+ii));
757	for ( jj = 1; jj<=ii+1; jj++ )
758	{
759	P=P+ran[1,jj]*m[jj];
760	}
761	chcoord[n-d+2+ii]=P;
762	P=0;
763	}
764	phi=r1,chcoord;
765	dbprint(printlevel-voice+2,"// (1 random change of coord.)");
766	I=simplify(lead(std(phi(i))),1);
767	attrib(I,"isSB",1);
768	nesttest=is_nested(I);
769	if ( nesttest == 1 )
770	{
771	break;
772	}
773	}
774	if ( nesttest == 0 )
775	{
776	"// WARNING from proc depthIdeal from lib mregular.lib:
777	// The procedure has entered in 30 loops and could not find a monomial
778	// ideal of nested type with the same depth as your ideal: in your
779	// example the method using random changes of coordinates may enter an
780	// infinite loop when the field is finite.
781	// Try removing this optional argument.";
782	return (-1);
783	}
784	}
785	}
786	//
787	// At this stage, we have obtained a monomial ideal I of nested type
788	// such that depth(S/i)=depth(S/I). We now compute depth(I).
789	//
790	//----- When S/i is Cohen-Macaulay:
791	for ( ii = n-d+2; ii <= n+1; ii++ )
792	{
793	K=K+select(I,ii);
794	}
795	if ( size(K) == 0 )
796	{
797	time=rtimer-time;
798	// Additional information:
799	dbprint(printlevel-voice+2,
800	"// Dimension of S/i: "+string(d)+" (S/i is Cohen-Macaulay)");
801	dbprint(printlevel-voice+2,
802	"// Time for computing depth: " + string(time) + " sec.");
803	dbprint(printlevel-voice+2,
804	"// The depth of S/i is:");
805	return(d);
806	}
807	//----- When d=1 (and S/i is not Cohen-Macaulay) ==> depth =0:
808	if ( d == 1 )
809	{
810	time=rtimer-time;
811	// Additional information:
812	dbprint(printlevel-voice+2,
813	"// Dimension of S/i: 1");
814	dbprint(printlevel-voice+2,
815	"// Time for computing depth: "+ string(time) + " sec.");
816	dbprint(printlevel-voice+2,
817	"// The depth of S/i is:");
818	return(0);
819
820	}
821	//----- Now d>1 and S/i is not Cohen-Macaulay:
822	//
823	//----- First, determine the last variable really occuring
824	lastv=n-d;
825	h=n;
826	while ( lastv == n-d and h > n-d )
827	{
828	K=select(I,h+1);
829	if ( size(K) == 0 )
830	{
831	h=h-1;
832	}
833	else
834	{
835	lastv=h;
836	}
837	}
838	//----- and compute the depth:
839	time=rtimer-time;
840	// Additional information:
841	dbprint(printlevel-voice+2,
842	"// Dimension of S/i: "+string(d));
843	dbprint(printlevel-voice+2,
844	"// Time for computing depth: "+ string(time) + " sec.");
845	dbprint(printlevel-voice+2,
846	"// The depth of S/i is:");
847	return(n-lastv);
848	}
849	example
850	{ "EXAMPLE:"; echo = 2;
851	ring r=0,(x,y,z,t,w),dp;
852	ideal i=y2t,x2y-x2z+yt2,x2y2,xyztw,x3z2,y5+xz3w-x2zw2,x7-yt2w4;
853	depthIdeal(i);
854	depthIdeal(lead(std(i)));
855	// Additional information is displayed if you change printlevel (=1);
856	}
857	////////////////////////////////////////////////////////////////////////////////
858	/*
859	Out-commented examples:
860	ring s=0,x(0..5),dp;
861	ideal i=x(2)^2-x(4)x(5),x(1)x(2)-x(0)x(5),x(0)x(2)-x(1)*x(4),
862	x(1)^2-x(3)x(5),x(0)x(1)-x(2)x(3),x(0)^2-x(3)x(4);
863	depthIdeal(i);
864	// Our procedure works when a min. graded free resol. can
865	// not be computed. In this easy example, depth can also
866	// be obtained using a m.g.f.r. (Auslander-Buchsbaum formula):
867	nvars(s)-ncols(betti(mres(i,0)))+1;
868	ring r1=0,(x,y,z,t),dp;
869	// Ex.2.5 in [Bermejo-Gimenez], Proc.Amer.Math.Soc. 128(5):
870	ideal i = x17y14-y31, x20y13, x60-y36z24-x20z20t20;
871	depthIdeal(i);
872	// Ex.2.9 in [Bermejo-Gimenez], Proc.Amer.Math.Soc. 128(5):
873	int k=43;
874	ideal j=x17y14-y31,x20y13,x60-y36z24-x20z20t20,y41z^k-y40z^(k+1);
875	depthIdeal(j);
876	// Example in Rk.2.10 in [Bermejo-Gimenez], ProcAMS 128(5):
877	ideal h=x2-3xy+5xt,xy-3y2+5yt,xz-3yz,2xt-yt,y2-yz-2yt;
878	depthIdeal(h);
879	// The initial ideal is not saturated
880	depthIdeal(lead(std(h)));
881	// More examples:
882	i=y4-t3z, x3t-y2z2, x3y2-t2z3, x6-tz5;
883	depthIdeal(i);
884	//
885	depthIdeal(maxideal(4));
886	//
887	ring r2=0,(x,y,z,t,w),dp;
888	ideal i = xy-zw,x3-yw2,x2z-y2w,y3-xz2,-y2z3+xw4+tw4+w5,-yz4+x2w3+xtw3+xw4,
889	-z5+x2tw2+x2w3+yw4;
890	depthIdeal(i);
891	//
892	ring r3=0,(x,y,z,t,w,u),dp;
893	ideal i=imap(r2,i);
894	depthIdeal(i);
895	// Next example is the defining ideal of the 2nd. Veronesean of P3, a variety
896	// in P8 which is arithmetically Cohen-Macaulay:
897	ring r4=0,(a,b,c,d,x(0..9)),dp;
898	ideal i= x(0)-ab,x(1)-ac,x(2)-ad,x(3)-bc,x(4)-bd,x(5)-cd,
899	x(6)-a2,x(7)-b2,x(8)-c2,x(9)-d2;
900	ideal ei=eliminate(i,abcd);
901	ring r5=0,x(0..9),dp;
902	ideal i=imap(r4,ei);
903	depthIdeal(i);
904	// Here is an example where the computation of a m.g.f.r. of I costs:
905	ring r8=0,(x,y,z,t,u,a,b),dp;
906	ideal i=u-b40,t-a40,x-a23b17,y-a22b18+ab39,z-a25b15;
907	ideal ei=eliminate(i,ab); // It takes a few seconds to compute the ideal
908	ring r9=0,(x,y,z,t,u),dp;
909	ideal i=imap(r8,ei);
910	depthIdeal(i); // This is very fast.
911	// Now you can use mres(i,0) to compute a m.g.f.r. of the ideal!
912	//
913	// Another one:
914	ring r10=0,(x(0..8),s,t),dp;
915	ideal i=x(0)-st24,x(1)-s2t23,x(2)-s3t22,x(3)-s9t16,x(4)-s11t14,x(5)-s18t7,
916	x(6)-s24t,x(7)-t25,x(8)-s25;
917	ideal ei=eliminate(i,st);
918	ring r11=0,x(0..8),dp;
919	ideal i=imap(r10,ei);
920	depthIdeal(i);
921	// More examples where not even sres works:
922	// Be careful: elimination takes some time here, but it succeeds!
923	ring r12=0,(s,t,u,x(0..14)),dp;
924	ideal i=x(0)-st6u8,x(1)-s5t3u7,x(2)-t11u4,x(3)-s9t4u2,x(4)-s2t7u6,x(5)-s7t7u,
925	x(6)-s10t5,x(7)-s4t6u5,x(8)-s13tu,x(9)-s14u,x(10)-st2u12,x(11)-s3t9u3,
926	x(12)-s15,x(13)-t15,x(14)-u15;
927	ideal ei=eliminate(i,stu);
928	size(ei);
929	ring r13=0,x(0..14),dp;
930	ideal i=imap(r12,ei);
931	size(i);
932	depthIdeal(i);
933	//
934	*/
935	///////////////////////////////////////////////////////////////////////////////
936	///////////////////////////////////////////////////////////////////////////////
937	///////////////////////////////////////////////////////////////////////////////
938
939	proc satiety (ideal i, list #)
940	"
941	USAGE: satiety (i[,e]); i ideal, e integer
942	RETURN: an integer, the satiety of i.
943	(returns -1 if i is not homogeneous)
944	ASSUME: i is a homogeneous ideal of the basering S=K[x(0)..x(n)].
945	e=0: (default)
946	The satiety is computed determining the fresh elements in the
947	socle of i. It works over arbitrary fields.
948	e=1: Makes random changes of coordinates to find a monomial ideal
949	with same satiety. It works over infinite fields only. If K
950	is finite, it works if it terminates, but may result in an
951	infinite loop. After 30 loops, a warning message is displayed
952	and -1 is returned.
953	THEORY: The satiety, or saturation index, of a homogeneous ideal i is the
954	least integer s such that, for all d>=s, the degree d part of the
955	ideals i and isat=sat(i,maxideal(1))[1] coincide.
956	NOTE: If printlevel > 0 (default = 0), dim(S/i) is also displayed.
957	EXAMPLE: example satiety; shows some examples
958	"
959	{
960	//--------------------------- initialisation ---------------------------------
961	int e,ii,jj,h,d,time,lastv,nesttest,NPtest,nl,sat;
962	intmat ran;
963	def r0 = basering;
964	int n = nvars(r0)-1;
965	if ( size(#) > 0 )
966	{
967	e = #[1];
968	}
969	string s = "ring r1 = ",charstr(r0),",x(0..n),dp;";
970	execute(s);
971	ideal i,sbi,I,K,chcoord,m,KK;
972	poly P;
973	map phi;
974	i = fetch(r0,i);
975	time=rtimer;
976	sbi=std(i);
977	//----- Check ideal homogeneous
978	if ( homog(sbi) == 0 )
979	{
980	"// WARNING from proc satiety from lib mregular.lib:
981	// The ideal is not homogeneous!";
982	return (-1);
983	}
984	I=simplify(lead(sbi),1);
985	attrib(I,"isSB",1);
986	d=dim(I);
987	//----- If the ideal i is not proper:
988	if ( d == -1 )
989	{
990	dbprint(printlevel-voice+2,
991	"// The ideal i is (1)!
992	// Its satiety is:");
993	return (0);
994	}
995	//----- If the ideal i is 0:
996	if ( size(I) == 0 )
997	{
998	dbprint(printlevel-voice+2,
999	"// The ideal i is (0)!
1000	// Its satiety is:");
1001	return (0);
1002	}
1003	//----- When the ideal i is 0-dimensional:
1004	if ( d == 0 )
1005	{
1006	sat=maxdeg1(minbase(quotient(I,maxideal(1))))+1;
1007	time=rtimer-time;
1008	// Additional information:
1009	dbprint(printlevel-voice+2,
1010	"// Dimension of S/i: 0");
1011	dbprint(printlevel-voice+2,
1012	"// Time for computing the satiety: " + string(time) + " sec.");
1013	dbprint(printlevel-voice+2,
1014	"// The satiety of i is:");
1015	return (sat);
1016	}
1017	//----- When one has option e=1:
1018	//
1019	//----- Determine the situation: NT, or NP, or nothing.
1020	//----- Choose the method depending on the situation in order to
1021	//----- get the mon. ideal of nested type associated to i
1022	if ( e == 1 )
1023	{
1024	NPtest=is_NP(I);
1025	if ( NPtest == 0 )
1026	{
1027	while ( nl < 30 )
1028	{
1029	chcoord=select1(maxideal(1),1..(n-d+1));
1030	nl=nl+1;
1031	for ( ii = 1; ii<=d; ii++ )
1032	{
1033	ran=random(100,1,n-d+ii);
1034	ran=intmat(ran,1,n-d+ii+1);
1035	ran[1,n-d+ii+1]=1;
1036	m=select1(maxideal(1),1..(n-d+1+ii));
1037	for ( jj = 1; jj<=n-d+ii+1; jj++ )
1038	{
1039	P=P+ran[1,jj]*m[jj];
1040	}
1041	chcoord[n-d+1+ii]=P;
1042	P=0;
1043	}
1044	phi=r1,chcoord;
1045	dbprint(printlevel-voice+2,"// (1 random change of coord.)");
1046	I=simplify(lead(std(phi(i))),1);
1047	attrib(I,"isSB",1);
1048	NPtest=is_NP(I);
1049	if ( NPtest == 1 )
1050	{
1051	break;
1052	}
1053	}
1054	if ( NPtest == 0 )
1055	{
1056	"// WARNING from proc satiety from lib mregular.lib:
1057	// The procedure has entered in 30 loops and could not put the variables
1058	// in Noether position: in your example the method using random changes
1059	// of coordinates may enter an infinite loop when the field is finite.
1060	// Try removing the optional argument.";
1061	return (-1);
1062	}
1063	i=phi(i);
1064	}
1065	nesttest=is_nested(I);
1066	if ( nesttest == 0 )
1067	{
1068	while ( nl < 30 )
1069	{
1070	chcoord=select1(maxideal(1),1..(n-d+2));
1071	nl=nl+1;
1072	for ( ii = 1; ii<=d-1; ii++ )
1073	{
1074	ran=random(100,1,ii);
1075	ran=intmat(ran,1,ii+1);
1076	ran[1,ii+1]=1;
1077	m=select1(maxideal(1),(n-d+2)..(n-d+2+ii));
1078	for ( jj = 1; jj<=ii+1; jj++ )
1079	{
1080	P=P+ran[1,jj]*m[jj];
1081	}
1082	chcoord[n-d+2+ii]=P;
1083	P=0;
1084	}
1085	phi=r1,chcoord;
1086	dbprint(printlevel-voice+2,"// (1 random change of coord.)");
1087	I=simplify(lead(std(phi(i))),1);
1088	attrib(I,"isSB",1);
1089	nesttest=is_nested(I);
1090	if ( nesttest == 1 )
1091	{
1092	break;
1093	}
1094	}
1095	if ( nesttest == 0 )
1096	{
1097	"// WARNING from proc satiety from lib mregular.lib:
1098	// The procedure has entered in 30 loops and could not find a monomial
1099	// ideal of nested type with the same satiety as your ideal: in your
1100	// example the method using random changes of coordinates may enter an
1101	// infinite loop when the field is finite.
1102	// Try removing the optional argument.";
1103	return (-1);
1104	}
1105	}
1106	//
1107	// At this stage, we have obtained a monomial ideal I of nested type
1108	// such that depth(S/i)=depth(S/I). We now compute depth(I).
1109	//
1110	//----- When S/i is Cohen-Macaulay:
1111	//
1112	for ( ii = n-d+2; ii <= n+1; ii++ )
1113	{
1114	K=K+select(I,ii);
1115	}
1116	if ( size(K) == 0 )
1117	{
1118	time=rtimer-time;
1119	// Additional information:
1120	dbprint(printlevel-voice+2,
1121	"// Dimension of S/i: "+string(d));
1122	dbprint(printlevel-voice+2,
1123	"// Time for computing satiety: " + string(time) + " sec.");
1124	dbprint(printlevel-voice+2,
1125	"// The satiety of i is:");
1126	return(0);
1127	}
1128	//----- When d=1 (and S/i is not Cohen-Macaulay) ==> depth =0:
1129	if ( d == 1 )
1130	{
1131	KK=simplify(reduce(quotient(I,maxideal(1)),I),2);
1132	sat=maxdeg1(KK)+1;
1133	time=rtimer-time;
1134	// Additional information:
1135	dbprint(printlevel-voice+2,
1136	"// Dimension of S/i: 1");
1137	dbprint(printlevel-voice+2,
1138	"// Time for computing satiety: "+ string(time) + " sec.");
1139	dbprint(printlevel-voice+2,
1140	"// The satiety of i is:");
1141	return(sat);
1142	}
1143	//----- Now d>1 and S/i is not Cohen-Macaulay:
1144	//
1145	//----- First, determine the last variable really occuring
1146	lastv=n-d;
1147	h=n;
1148	while ( lastv == n-d and h > n-d )
1149	{
1150	K=select(I,h+1);
1151	if ( size(K) == 0 )
1152	{
1153	h=h-1;
1154	}
1155	else
1156	{
1157	lastv=h;
1158	}
1159	}
1160	//----- and compute the satiety:
1161	sat=0;
1162	if ( lastv == n )
1163	{
1164	KK=simplify(reduce(quotient(I,maxideal(1)),I),2);
1165	sat=maxdeg1(KK)+1;
1166	}
1167	time=rtimer-time;
1168	// Additional information:
1169	dbprint(printlevel-voice+2,
1170	"// Dimension of S/i: "+string(d));
1171	dbprint(printlevel-voice+2,
1172	"// Time for computing satiety: "+ string(time) + " sec.");
1173	dbprint(printlevel-voice+2,
1174	"// The satiety of i is:");
1175	return(sat);
1176	}
1177	//---- If no option: direct computation
1178	sat=maxdeg1(reduce(quotient(i,maxideal(1)),sbi))+1;
1179	time=rtimer-time;
1180	// Additional information:
1181	dbprint(printlevel-voice+2,
1182	"// Dimension of S/i: "+string(d)+";");
1183	dbprint(printlevel-voice+2,
1184	"// Time for computing satiety: "+ string(time) + " sec.");
1185	dbprint(printlevel-voice+2,
1186	"// The satiety of i is:");
1187	return(sat);
1188	}
1189	example
1190	{ "EXAMPLE:"; echo = 2;
1191	ring r=0,(x,y,z,t,w),dp;
1192	ideal i=y2t,x2y-x2z+yt2,x2y2,xyztw,x3z2,y5+xz3w-x2zw2,x7-yt2w4;
1193	satiety(i);
1194	ideal I=lead(std(i));
1195	satiety(I); // First method: direct computation
1196	satiety(I,1); // Second method: doing changes of coordinates
1197	// Additional information is displayed if you change printlevel (=1);
1198	}
1199	////////////////////////////////////////////////////////////////////////////////
1200	/*
1201	Out-commented examples:
1202	ring s1=0,(x,y,z,t),dp;
1203	ideal I=zt3,z2t2,yz2t,xz2t,xy2t,x3y;
1204	satiety(I);
1205	satiety(I,1);
1206	// Another example:
1207	ring s2=0,(z,y,x),dp;
1208	ideal I=z38,z26y2,z14y4,z12x,z10x5,z8x9,z6x16,z4x23,z2y6,y32;
1209	satiety(I);
1210	satiety(I,1);
1211	// One more:
1212	ring s3=0,(s,t,u,x(0..8)),dp;
1213	ideal i=x(0)-st6u8,x(1)-s5t3u7,x(2)-t11u4,x(3)-s9t4u2,
1214	x(4)-s2t7u6,x(5)-s7t7u,x(6)-s15,x(7)-t15,x(8)-u15;
1215	ideal ei=eliminate(i,stu);
1216	size(ei);
1217	ring s4=0,x(0..8),dp;
1218	ideal i=imap(s3,ei);
1219	ideal m=maxideal(1);
1220	m[8]=m[8]+m[7];
1221	map phi=m;
1222	ideal phii=phi(i);
1223	ideal nI=lead(std(phii));
1224	ring s5=0,x(0..7),dp;
1225	ideal nI=imap(s4,nI);
1226	satiety(nI);
1227	satiety(nI,1);
1228	ideal I1=subst(nI,x(7),1);
1229	ring s6=0,x(0..6),dp;
1230	ideal I1=imap(s5,I1);
1231	satiety(I1);
1232	satiety(I1,1);
1233	ideal I2=subst(I1,x(6),1);
1234	ring s7=0,x(0..5),dp;
1235	ideal I2=imap(s6,I2);
1236	satiety(I2);
1237	satiety(I2,1);
1238	//
1239	*/
1240	///////////////////////////////////////////////////////////////////////////////
1241	///////////////////////////////////////////////////////////////////////////////
1242	///////////////////////////////////////////////////////////////////////////////
1243
1244	proc regMonCurve (list #)
1245	"
1246	USAGE: regMonCurve (a0,...,an) ; ai integers with a0=0 < a1 < ... < an=:d
1247	RETURN: an integer, the Castelnuovo-Mumford regularity of the projective
1248	monomial curve C in Pn(K) parametrically defined by
1249	x(0) = t^d , x(1) = s^(a1)t^(d-a1) , ..... , x(n) = s^d
1250	where K is the field of complex numbers.
1251	(returns -1 if a0=0 < a1 < ... < an is not satisfied)
1252	ASSUME: a0=0 < a1 < ... < an are integers.
1253	NOTES: 1. The defining ideal of the curve C, I in S=K[x(0),...,x(n)], is
1254	determined by elimination.
1255	2. The procedure regIdeal has been improved in this case since one
1256	knows beforehand that the monomial ideal J=lead(std(I)) is of
1257	nested type if the monomial ordering is dp, and that
1258	reg(C)=reg(J) (see preprint 'Saturation and Castelnuovo-Mumford
1259	regularity' by Bermejo-Gimenez, 2004).
1260	3. If printlevel > 0 (default = 0) additional info is displayed:
1261	- It says whether C is arithmetically Cohen-Macaulay or not.
1262	- If C is not arith. Cohen-Macaulay, end(H^1(S/I)) is computed
1263	and an upper bound for the a-invariant of S/I is given.
1264	- It also determines one step of the minimal graded free
1265	resolution (m.g.f.r.) of I where the regularity is attained
1266	and gives the value of the regularity of the Hilbert function
1267	of S/I when reg(I) is attained at the last step of a m.g.f.r.
1268	EXAMPLE: example regMonCurve; shows some examples
1269	"
1270	{
1271	//--------------------------- initialisation ---------------------------------
1272	int ii,H,h,hh,time,ttime,firstind,lastind;
1273	int n = size(#)-1;
1274	//------------------ Check assumptions on integers -------------------------
1275	if ( #[1] != 0 )
1276	{"// WARNING from proc regMonCurve from lib mregular.lib:
1277	// USAGE: your input must be a list of integers a0,a1,...,an such that
1278	// a0=0 < a1 < a2 < ... < an";
1279	return(-1);
1280	}
1281	for ( ii=1; ii<= n; ii++ )
1282	{
1283	if ( #[ii] >= #[ii+1] )
1284	{
1285	"// WARNING from proc regMonCurve from lib mregular.lib:
1286	// USAGE: your input must be a list of integers a0,a1,...,an such that
1287	// a0=0 < a1 < a2 < ... < an";
1288	return(-1);
1289	}
1290	}
1291	ring R=0,(x(0..n),s,t),dp;
1292	ideal param,m,i;
1293	poly f(0..n);
1294	for (ii=0;ii<=n;ii++)
1295	{
1296	f(ii)=s^(#[n+1]-#[ii+1])*t^(#[ii+1]);
1297	param=param+f(ii);
1298	}
1299	m=subst(maxideal(1),s,0);
1300	m=simplify(subst(m,t,0),2);
1301	i=matrix(m)-matrix(param);
1302	ttime=rtimer;
1303	i=eliminate(i,st);
1304	ring r=0,(x(1..n),x(0)),dp;
1305	ideal i,I;
1306	i=imap(R,i);
1307	I=minbase(lead(std(i)));
1308	attrib(I,"isSB",1);
1309	ttime=rtimer-ttime;
1310	time=rtimer;
1311	ring nr=0,x(1..n),dp;
1312	ideal I,K,KK,J;
1313	I=imap(r,I);
1314	attrib(I,"isSB",1);
1315	K=select(I,n);
1316	//------------------ Cohen-Macaulay case ------------
1317	if ( size(K) == 0 )
1318	{
1319	ring mr=0,x(1..n-1),dp;
1320	ideal I=imap(nr,I);
1321	H=maxdeg1(minbase(quotient(I,maxideal(1))))+1;
1322	time=rtimer-time;
1323	// Additional information:
1324	dbprint(printlevel-voice+2,
1325	"// The sequence of integers defines a monomial curve C in P"
1326	+ string(n));
1327	dbprint(printlevel-voice+2,
1328	"// C is arithmetically Cohen-Macaulay");
1329	dbprint(printlevel-voice+2,
1330	"// Regularity attained at the last step of a m.g.f.r. of I(C)");
1331	dbprint(printlevel-voice+2,
1332	"// Regularity of the Hilbert function of S/I(C): "
1333	+ string(H-2));
1334	dbprint(printlevel-voice+2,
1335	"// Time for computing ideal I(C) (by elimination): "
1336	+ string(ttime) + " sec.");
1337	dbprint(printlevel-voice+2,
1338	"// Time for computing reg(C) once I(C) has been determined: "
1339	+ string(time) + " sec.");
1340	dbprint(printlevel-voice+2,
1341	"// The Castelnuovo-Mumford regularity of C is:");
1342	return(H);
1343	}
1344	else
1345	{
1346	KK=simplify(reduce(quotient(I,maxideal(1)),I),2);
1347	firstind=maxdeg1(KK)+1;
1348	J=subst(I,x(n),1);
1349	ring mr=0,x(1..n-1),dp;
1350	ideal J=imap(nr,J);
1351	lastind=maxdeg1(minbase(quotient(J,maxideal(1))))+1;
1352	H=firstind;
1353	if ( lastind > H )
1354	{
1355	H=lastind;
1356	}
1357	time=rtimer-time;
1358	// Additional information:
1359	dbprint(printlevel-voice+2,
1360	"// The sequence of integers defines a monomial curve C in P"
1361	+ string(n));
1362	dbprint(printlevel-voice+2,
1363	"// C is not arithmetically Cohen-Macaulay");
1364	dbprint(printlevel-voice+2,
1365	"// end(H^1(S/I(C))) = "
1366	+string(firstind-2));
1367	dbprint(printlevel-voice+2,
1368	"// Upper bound for the a-invariant of S/I(C): end(H^2(S/I(C))) <= "
1369	+string(lastind-3));
1370	if ( H == firstind )
1371	{
1372	dbprint(printlevel-voice+2,
1373	"// Regularity attained at the last step of a m.g.f.r. of I(C)");
1374	dbprint(printlevel-voice+2,
1375	"// Regularity of the Hilbert function of S/I(C): "
1376	+ string(H-1));
1377	}
1378	else
1379	{
1380	dbprint(printlevel-voice+2,
1381	"// Regularity attained at the second last step of a m.g.f.r. of I(C)");
1382	dbprint(printlevel-voice+2,
1383	"// (and not attained at the last step)");
1384	}
1385	dbprint(printlevel-voice+2,
1386	"// Time for computing ideal I(C) (by elimination): "
1387	+ string(ttime) + " sec.");
1388	dbprint(printlevel-voice+2,
1389	"// Time for computing reg(C) once I(C) has been determined: "
1390	+ string(time) + " sec.");
1391	dbprint(printlevel-voice+2,
1392	"// The Castelnuovo-Mumford regularity of C is:");
1393	return(H);
1394	}
1395	}
1396	example
1397	{ "EXAMPLE:"; echo = 2;
1398	// The 1st example is the twisted cubic:
1399	regMonCurve(0,1,2,3);
1400	// The 2nd. example is the non arithm. Cohen-Macaulay monomial curve in P4
1401	// parametrized by: x(0)-s6,x(1)-s5t,x(2)-s3t3,x(3)-st5,x(4)-t6:
1402	regMonCurve(0,1,3,5,6);
1403	// Additional information is displayed if you change printlevel (=1);
1404	}
1405	////////////////////////////////////////////////////////////////////////////////
1406	/*
1407	Out-commented examples:
1408	//
1409	// The sequence of integers must be strictly increasing
1410	regMonCurve(1,4,6,9);
1411	regMonCurve(0,3,8,5,23);
1412	regMonCurve(0,4,7,7,9);
1413	//
1414	// A curve in P3 s.t. the regularity is attained at the last step:
1415	regMonCurve(0,2,12,15);
1416	//
1417	// A curve in P4 s.t. the regularity attained at the last but one
1418	// but NOT at the last step (Ex. 3.3 Preprint 2004):
1419	regMonCurve(0,5,9,11,20);
1420	//
1421	// A curve in P8 s.t. the m.g.f.r. of the defining ideal is not easily
1422	// obtained through m.g.f.r.:
1423	regMonCurve(0,1,2,3,9,11,18,24,25);
1424	//
1425	// A curve in P11 of degree 37:
1426	regMonCurve(0,1,2,7,16,17,25,27,28,30,36,37);
1427	// It takes some time to compute the eliminated ideal; the computation of
1428	// the regularity is then rather fast as one can check using proc regIdeal:
1429	ring q=0,(s,t,x(0..11)),dp;
1430	ideal i=x(0)-st36,x(1)-s2t35,x(2)-s7t30,x(3)-s16t21,x(4)-s17t20,x(5)-s25t12,
1431	x(6)-s27t10,x(7)-s28t9,x(8)-s30t7,x(9)-s36t,x(10)-s37,x(11)-t37;
1432	ideal ei=eliminate(i,st);
1433	ring qq=0,x(0..11),dp;
1434	ideal i=imap(q,ei);
1435	regIdeal(i);
1436	//
1437	// A curve in P14 of degree 55:
1438	regMonCurve(0,1,2,7,16,17,25,27,28,30,36,37,40,53,55);
1439	//
1440	*/
1441	///////////////////////////////////////////////////////////////////////////////
1442	///////////////////////////////////////////////////////////////////////////////
1443	///////////////////////////////////////////////////////////////////////////////
1444
1445	proc NoetherPosition (ideal i)
1446	"
1447	USAGE: NoetherPosition (i); i ideal
1448	RETURN: ideal such that, for the homogeneous linear transformation
1449	map phi=S,NoetherPosition(i);
1450	one has that K[x(n-d+1),...,x(n)] is a Noether normalization of
1451	S/phi(i) where S=K[x(0),...x(n)] is the basering and d=dim(S/i).
1452	(returns -1 if i = (0) or (1)).
1453	ASSUME: The field K is infinite and i is a nonzero proper ideal.
1454	NOTE: 1. It works also if K is a finite field if it terminates, but
1455	may result in an infinite loop. If the procedure enters more
1456	than 30 loops, -1 is returned and a warning message is displayed.@*
1457	2. If printlevel > 0 (default = 0), additional info is displayed:
1458	dim(S/i) and K[x(n-d+1),...,x(n)] are given.
1459	EXAMPLE: example NoetherPosition; shows some examples
1460	"
1461	{
1462	//--------------------------- initialisation ---------------------------------
1463	int ii,jj,d,time,nl,NPtest;
1464	intmat ran;
1465	def r0 = basering;
1466	ideal K,chcoord;
1467	int n = nvars(r0)-1;
1468	string s = "ring r1 = ",charstr(r0),",x(0..n),dp;";
1469	execute(s);
1470	ideal i,sbi,I,K,chcoord,m;
1471	poly P;
1472	map phi;
1473	i = fetch(r0,i);
1474	time=rtimer;
1475	sbi=std(i);
1476	I=simplify(lead(sbi),1);
1477	attrib(I,"isSB",1);
1478	d=dim(I);
1479	//----- If the ideal i is not proper:
1480	if ( d == -1 )
1481	{
1482	"// WARNING from proc NoetherPosition from lib mregular.lib:
1483	// The ideal i is (1)!";
1484	return (-1);
1485	}
1486	//----- If the ideal i is 0:
1487	if ( size(I) == 0 )
1488	{
1489	"// WARNING from proc NoetherPosition from lib mregular.lib:
1490	// The ideal i is (0)!";
1491	return (-1);
1492	}
1493	//----- When the ideal i is 0-dimensional:
1494	if ( d == 0 )
1495	{
1496	time=rtimer-time;
1497	// Additional information:
1498	dbprint(printlevel-voice+2,
1499	"// Dimension of S/i: 0");
1500	dbprint(printlevel-voice+2,
1501	"// Time for computing a Noether normalization: "
1502	+ string(time) + " sec.");
1503	dbprint(printlevel-voice+2,
1504	"// K is a Noether normalization of S/phi(i)");
1505	dbprint(printlevel-voice+2,
1506	"// where the map phi: S --> S is:");
1507	setring r0;
1508	return (maxideal(1));
1509	}
1510	NPtest=is_NP(I);
1511	if ( NPtest == 1 )
1512	{
1513	K=x(n-d+1..n);
1514	setring r0;
1515	K=fetch(r1,K);
1516	time=rtimer-time;
1517	// Additional information:
1518	dbprint(printlevel-voice+2,
1519	"// Dimension of S/i: " + string(d) );
1520	dbprint(printlevel-voice+2,
1521	"// Time for computing a Noether normalization: " +
1522	string(time) + " sec.");
1523	dbprint(printlevel-voice+2,
1524	"// K[" + string(K) +
1525	"] is a Noether normalization of S/phi(i)");
1526	dbprint(printlevel-voice+2,
1527	"// where the map phi: S --> S is:");
1528	return (maxideal(1));
1529	}
1530	//---- Otherwise, random change of coordinates and
1531	//---- test for Noether normalization.
1532	//---- If we were unlucky, another change of coord. will be done:
1533	while ( nl < 30 )
1534	{
1535	chcoord=select1(maxideal(1),1..(n-d+1));
1536	nl=nl+1;
1537	for ( ii = 1; ii<=d; ii++ )
1538	{
1539	ran=random(100,1,n-d+ii);
1540	ran=intmat(ran,1,n-d+ii+1);
1541	ran[1,n-d+ii+1]=1;
1542	m=select1(maxideal(1),1..(n-d+1+ii));
1543	for ( jj = 1; jj<=n-d+ii+1; jj++ )
1544	{
1545	P=P+ran[1,jj]*m[jj];
1546	}
1547	chcoord[n-d+1+ii]=P;
1548	P=0;
1549	}
1550	phi=r1,chcoord;
1551	dbprint(printlevel-voice+2,"// (1 random change of coord.)");
1552	I=simplify(lead(std(phi(i))),1);
1553	attrib(I,"isSB",1);
1554	NPtest=is_NP(I);
1555	if ( NPtest == 1 )
1556	{
1557	K=x(n-d+1..n);
1558	setring r0;
1559	K=fetch(r1,K);
1560	chcoord=fetch(r1,chcoord);
1561	time=rtimer-time;
1562	// Additional information:
1563	dbprint(printlevel-voice+2,
1564	"// Dimension of S/i: " + string(d) );
1565	dbprint(printlevel-voice+2,
1566	"// Time for computing a Noether normalization: " +
1567	string(time) + " sec.");
1568	dbprint(printlevel-voice+2,
1569	"// K[" + string(K) +
1570	"] is a Noether normalization of S/phi(i)");
1571	dbprint(printlevel-voice+2,
1572	"// where the map phi: S --> S is:");
1573	return (chcoord);
1574	}
1575	}
1576	"// WARNING from proc NoetherPosition from lib mregular.lib:
1577	// The procedure has entered in more than 30 loops: in your example
1578	// the method may enter an infinite loop over a finite field!";
1579	return (-1);
1580	}
1581	example
1582	{ "EXAMPLE:"; echo = 2;
1583	ring r=0,(x,y,z,t,u),dp;
1584	ideal i1=y,z,t,u; ideal i2=x,z,t,u; ideal i3=x,y,t,u; ideal i4=x,y,z,u;
1585	ideal i5=x,y,z,t; ideal i=intersect(i1,i2,i3,i4,i5);
1586	map phi=r,NoetherPosition(i);
1587	phi;
1588	ring r5=5,(x,y,z,t,u),dp;
1589	ideal i=imap(r,i);
1590	map phi=r5,NoetherPosition(i);
1591	phi;
1592	// Additional information is displayed if you change printlevel (=1);
1593	}
1594	///////////////////////////////////////////////////////////////////////////////
1595	///////////////////////////////////////////////////////////////////////////////
1596
1597	proc is_NP (ideal i)
1598	"
1599	USAGE: is_NP (i); i ideal
1600	RETURN: 1 if K[x(n-d+1),...,x(n)] is a Noether normalization of
1601	S/i where S=K[x(0),...x(n)] is the basering, and d=dim(S/i),
1602	0 otherwise.
1603	(returns -1 if i=(0) or i=(1)).
1604	ASSUME: i is a nonzero proper homogeneous ideal.
1605	NOTE: 1. If i is not homogeneous and is_NP(i)=1 then K[x(n-d+1),...,x(n)]
1606	is a Noether normalization of S/i. The converse may be wrong if
1607	the ideal is not homogeneous.
1608	2. is_NP is used in the procedures regIdeal, depthIdeal, satiety,
1609	and NoetherPosition.
1610	EXAMPLE: example is_NP; shows some examples
1611	"
1612	{
1613	//--------------------------- initialisation ---------------------------------
1614	int ii,d,dz;
1615	def r0 = basering;
1616	int n = nvars(r0)-1;
1617	string s = "ring r1 = ",charstr(r0),",x(0..n),dp;";
1618	execute(s);
1619	ideal i,sbi,I,J;
1620	i = fetch(r0,i);
1621	sbi=std(i);
1622	I=simplify(lead(sbi),1);
1623	attrib(I,"isSB",1);
1624	d=dim(I);
1625	//----- If the ideal i is not proper:
1626	if ( d == -1 )
1627	{
1628	"// WARNING from proc is_NP from lib mregular.lib:
1629	// The ideal i is (1)!";
1630	return (-1);
1631	}
1632	//----- If the ideal i is 0:
1633	if ( size(I) == 0 )
1634	{
1635	"// WARNING from proc is_NP from lib mregular.lib:
1636	// The ideal i is (0)!";
1637	return (-1);
1638	}
1639	//----- When the ideal i is 0-dimensional:
1640	if ( d == 0 )
1641	{
1642	return (1);
1643	}
1644	//----- Check Noether position
1645	J=I;
1646	for ( ii = n-d+1; ii <= n; ii++ )
1647	{
1648	J=subst(J,x(ii),0);
1649	}
1650	attrib(J,"isSB",1);
1651	dz=dim(J);
1652	if ( dz == d )
1653	{
1654	return (1);
1655	}
1656	else
1657	{
1658	return(0);
1659	}
1660	}
1661	example
1662	{ "EXAMPLE:"; echo = 2;
1663	ring r=0,(x,y,z,t,u),dp;
1664	ideal i1=y,z,t,u; ideal i2=x,z,t,u; ideal i3=x,y,t,u; ideal i4=x,y,z,u;
1665	ideal i5=x,y,z,t; ideal i=intersect(i1,i2,i3,i4,i5);
1666	is_NP(i);
1667	ideal ch=x,y,z,t,x+y+z+t+u;
1668	map phi=ch;
1669	is_NP(phi(i));
1670	}
1671	///////////////////////////////////////////////////////////////////////////////
1672	///////////////////////////////////////////////////////////////////////////////
1673
1674	proc is_nested (ideal i)
1675	"
1676	USAGE: is_nested (i); i monomial ideal
1677	RETURN: 1 if i is of nested type, 0 otherwise.
1678	(returns -1 if i=(0) or i=(1)).
1679	ASSUME: i is a nonzero proper monomial ideal.
1680	NOTES: 1. The ideal must be monomial, otherwise the result has no meaning
1681	(so check this before using this procedure).@*
1682	2. is_nested is used in procedures depthIdeal, regIdeal and satiety.@*
1683	3. When i is a monomial ideal of nested type of S=K[x(0)..x(n)],
1684	the a-invariant of S/i coincides with the upper bound obtained
1685	using the procedure regIdeal with printlevel > 0.
1686	THEORY: A monomial ideal is of nested type if its associated primes are all
1687	of the form (x(0),...,x(i)) for some i<=n.
1688	(see definition and effective criterion to check this property in
1689	the preprint 'Saturation and Castelnuovo-Mumford regularity' by
1690	Bermejo-Gimenez, 2004).
1691	EXAMPLE: example is_nested; shows some examples
1692	"
1693	{
1694	//--------------------------- initialisation ---------------------------------
1695	int ii,d,tev,lastv,h,NPtest;
1696	def r0 = basering;
1697	int n = nvars(r0)-1;
1698	string s = "ring r1 = ",charstr(r0),",x(0..n),dp;";
1699	execute(s);
1700	ideal I,K,KK,LL;
1701	I = fetch(r0,i);
1702	I=minbase(I);
1703	attrib(I,"isSB",1);
1704	d=dim(I);
1705	//----- If the ideal i is not proper:
1706	if ( d == -1 )
1707	{
1708	"// WARNING from proc is_nested from lib mregular.lib:
1709	// The ideal i is (1)!";
1710	return (-1);
1711	}
1712	//----- If the ideal i is 0:
1713	if ( size(I) == 0 )
1714	{
1715	"// WARNING from proc is_nested from lib mregular.lib:
1716	// The ideal i is (0)!";
1717	return (-1);
1718	}
1719	//----- When the ideal i is 0-dimensional:
1720	if ( d == 0 )
1721	{
1722	return (1);
1723	}
1724	//----- Check Noether position
1725	NPtest=is_NP(I);
1726	if ( NPtest != 1 )
1727	{
1728	return (0);
1729	}
1730	//----- When ideal is 1-dim. + var. in Noether position -> Nested Type
1731	if ( d == 1 )
1732	{
1733	return (1);
1734	}
1735	//----- Determ. of the last variable really occuring
1736	lastv=n-d;
1737	h=n;
1738	while ( lastv == n-d and h > n-d )
1739	{
1740	K=select(I,h+1);
1741	if ( size(K) == 0 )
1742	{
1743	h=h-1;
1744	}
1745	else
1746	{
1747	lastv=h;
1748	}
1749	}
1750	//----- Check the second property by evaluation when NP + d>1
1751	KK=subst(I,x(lastv),1);
1752	for ( ii = n-lastv; ii<=d-2; ii++ )
1753	{
1754	LL=minbase(subst(I,x(n-ii-1),1));
1755	attrib(LL,"isSB",1);
1756	tev=size(reduce(KK,LL));
1757	if ( tev > 0 )
1758	{
1759	return(0);
1760	}
1761	KK=LL;
1762	}
1763	return(1);
1764	}
1765	example
1766	{ "EXAMPLE:"; echo = 2;
1767	ring s=0,(x,y,z,t),dp;
1768	ideal i1=x2,y3; ideal i2=x3,y2,z2; ideal i3=x3,y2,t2;
1769	ideal i=intersect(i1,i2,i3);
1770	is_nested(i);
1771	ideal ch=x,y,z,z+t;
1772	map phi=ch;
1773	ideal I=lead(std(phi(i)));
1774	is_nested(I);
1775	}
1776	///////////////////////////////////////////////////////////////////////////////
1777	///////////////////////////////////////////////////////////////////////////////
1778

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

Download in other formats:

Original Format