Context Navigation

mregular.lib @ 0dd77c2

spielwiese

Last change on this file since 0dd77c2 was 867e1a3, checked in by Hans Schönemann <hannes@…>, 14 years ago
no - for ideal/module: track 114 git-svn-id: file:///usr/local/Singular/svn/trunk@12387 2c84dea3-7e68-4137-9b89-c4e89433aadc
Property mode set to `100644`
File size: 57.0 KB

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

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

Download in other formats:

Original Format