Context Navigation

mregular.lib @ e182c8

spielwiese

Last change on this file since e182c8 was e182c8, checked in by Hans Schönemann <hannes@…>, 19 years ago
*lossen/greuel/hannes: new libs git-svn-id: file:///usr/local/Singular/svn/trunk@7847 2c84dea3-7e68-4137-9b89-c4e89433aadc
Property mode set to `100644`
File size: 53.9 KB

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

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

Download in other formats:

Original Format