Context Navigation

binresol.lib @ 380a17b

spielwiese

Last change on this file since 380a17b was 380a17b, checked in by Hans Schoenemann <hannes@…>, 11 years ago
fix: new version numbers for libs
Property mode set to `100644`
File size: 86.0 KB

Line
1	///////////////////////////////////////////////////////////////////////////
2	version="version binresol.lib 4.0.0.0 Jun_2013 ";
3	category="Commutaive algebra";
4	info="
5	LIBRARY: binresol.lib Combinatorial algorithm of resolution of singularities
6	of binomial ideals in arbitrary characteristic.
7	Binomial resolution algorithm of Blanco
8
9	AUTHORS: R. Blanco, rblanco@modulor.arq.uva.es,
10	@* G. Pfister, pfister@mathematik.uni-kl.de
11
12	PROCEDURES:
13	Eresol(J); computes a E-resolution of singularities of (J) in char 0
14	BINresol(J); computes a E-resolution of singularities of (J) (THE SECOND PART IS NOT IMPLEMENTED YET)
15	determinecenter(L1,L2,c,n,Y,a,mb,flag,control3); computes the next blowing-up center
16	Blowupcenter(L1,id,m,L2,c,n,h); makes the blowing-up
17	Nonhyp(Coef,expJ,sJ,n,flag,sums); computes the ideal generated by the non hyperbolic generators of expJ
18	inidata(K,k); verifies input data, a binomial ideal K of k generators
19	identifyvar(); identifies status of variables
20	data(K,k,n); transforms data on lists of length n
21	Edatalist(Coef,Exp,k,n,flag); gives the E-order of each term in Exp
22	EOrdlist(Coef,Exp,k,n,flag); computes the E-order of an ideal (giving in the language of lists)
23	maxEord(Coef,Exp,k,n,flag); computes de maximum E-order of an ideal given by Coef and Exp
24	ECoef(Coef,expP,sP,V,auxc,n,flag); Computes a simplified version of the E-Coeff ideal.
25	elimrep(L); removes repeated terms from a list
26	Emaxcont(Coef,Exp,k,n,flag); computes a list of hypersurfaces of E-maximal contact
27	cleanunit(mon,n,flag); clean the units in a monomial mon
28	resfunction(t,auxinv,nchart,n); composes the E-resolution function
29	calculateI(Coef,J,c,n,Y,a,b,D); computes the order of the non monomial part of an ideal J
30	Maxord(L,n); computes the maximum exponent of an exceptional monomial ideal
31	Gamma(L,c,n); computes the Gamma function for an exceptional monomial ideal given by L
32	convertdata(C,L,n,flag); computes the ideal corresponding to C,L
33	tradblwup(blwhist,n,Y,a,num); composes the blowing up at this chart
34	lcmofall(nchart,mobile); computes the lcm of the denominators of the E-orders for all the charts
35	computemcm(Eolist); computes the lcm of the denominators of the E-orders for one chart
36	constructH(Hhist,n,flag); construct the list of exceptional divisors accumulated at this chart
37	constructblwup(blwhist,n,chy,flag); construct the ideal defining the composition map
38	constructlastblwup(blwhist,n,chy,flag); construct the ideal defining the last blowup leading to this chart
39	genoutput(chart,mobile,nchart,nsons,n,q); generates the output for visualization
40	salida(idchart,chart,mobile,numson,previousa,n,q); generates the output for one chart
41	iniD(n); creates a list of lists of zeros of size n
42	sumlist(L1,L2); sums two lists component to component
43	reslist(L1,L2); subtracts two lists component to component
44	multiplylist(L,a); multiplies a list by a number, component to component
45	dividelist(L1,L2); divides two lists component to component
46	createlist(L1,L2); creates a list of lists of two elements
47	list0(n); creates a list of zeros of size n
48	";
49
50	LIB "general.lib";
51	LIB "qhmoduli.lib";
52	LIB "inout.lib";
53	LIB "poly.lib";
54	LIB "resolve.lib";
55	LIB "reszeta.lib";
56	LIB "resgraph.lib";
57	////////////////////////////////////////////////////////////////////////////
58
59	proc inidata(ideal K,int k)
60	"USAGE: inidata(K,k); K any ideal, k integer (!=0)
61	COMPUTE: Verifies the input data
62	RETURN: flag indicating if the ideal is binomial or not
63	EXAMPLE: example inidata; shows an example
64	"
65	{
66	int i;
67	for (i=1;i<=k; i++)
68	{
69	if (size(K[i])>2){return(0);}
70	}
71	return(1);
72	}
73	example
74	{"EXAMPLE:"; echo = 2;
75	ring r = 0,(x(1..3)),dp;
76	ideal J1=x(1)^4*x(2)^2, x(1)^2+x(3)^3;
77	inidata(J1,2);
78
79	ideal J2=x(1)^4*x(2)^2, x(1)^2+x(2)^3+x(3)^5;
80	inidata(J2,2);
81	}
82	/////////////////////////////////////////////////////////////////////////////////
83
84	proc identifyvar()
85	"USAGE: identifyvar();
86	COMPUTE: Asign 0 to variables x and 1 to variables y, only necessary at the beginning
87	RETURN: list, say l, of size the dimension of the basering
88	l[i] is: 0 if the i-th variable is x(i),
89	1 if the i-th variable is y(i)
90	EXAMPLE: example identifyvar; shows an example
91	"
92	{
93	int i,n;
94	list flaglist;
95
96	n=nvars(basering);
97	flaglist=list0(n);
98
99	for (i=1;i<=n; i++){if (varstr(i)[1]=="y"){flaglist[i]=1;}}
100
101	return(flaglist);
102	}
103	example
104	{"EXAMPLE:"; echo = 2;
105	ring r = 0,(x(1),y(2),x(3),y(4),x(5..7),y(8)),dp;
106	identifyvar();
107	}
108	////////////////////////////////////////////////////////////////////////////////
109
110	proc data(ideal K,int k,int n)
111	"USAGE: data(K,k,n); K any ideal, k integer (!=0), n integer (!=0)
112	COMPUTE: Construcs a list with the coefficients and exponents of one ideal
113	RETURN: lists of coefficients and exponents of K
114	EXAMPLE: example data; shows an example
115	"
116	{int i,j,lon;
117	number aa;
118	intvec cc;
119	list bb,dd,aux,ddaux,Coef,Exp;
120
121	for (i=1;i<=k; i++)
122	{ lon=size(K[i]);
123
124	// binomial
125	if (lon==2){aa=leadcoef(K[i][1]);
126	bb=aa;
127	Coef[i]=bb; // coefficients
128	cc=leadexp(K[i][1]); // exponents
129
130	// cc is an intvec, transform cc in dd, a list of lists
131	dd=cc[1..n];
132	aux[1]=dd;
133	// the same for the second term
134
135	aa=leadcoef(K[i][2]);
136	bb=aa;
137	Coef[i]=Coef[i] + bb; // all the coefficients of i-th generator of K
138	cc=leadexp(K[i][2]);
139
140	dd=cc[1..n];
141	aux[2]=dd;
142	Exp[i]=aux;}
143
144	// monomial
145	if (lon==1){aux=list();
146	aa=leadcoef(K[i][1]);
147	bb=aa;
148	Coef[i]=bb;
149	cc=leadexp(K[i][1]);
150	dd=cc[1..n];
151	aux[1]=dd;
152	Exp[i]=aux;}
153	} //end for
154	return(Coef,Exp);
155	}
156	example
157	{"EXAMPLE:"; echo = 2;
158	ring r = 0,(x(1..3)),dp;
159	ideal J=x(1)^4*x(2)^2, x(1)^2-x(3)^3;
160	data(J,2,3);
161	}
162	//////////////////////////////////////////////////////
163
164	proc Edatalist(list Coef,list Exp,int k,int n,list flaglist)
165	"USAGE: Edatalist(Coef,Exp,k,n,flaglist);
166	Coef,Exp,flaglist lists, k,n, integers
167	Exp is a list of lists of exponents, k=size(Exp)
168	COMPUTE: computes a list with the E-order of each term
169	RETURN: a list with the E-order of each term
170	EXAMPLE: example Edatalist; shows an example
171	"
172	{int i,j,lon,mm;
173	list dd,ss,sums;
174	number aux,aux1,aux2;
175
176	for (i=1;i<=k;i++)
177	{
178	lon=size(Coef[i]);
179	if (lon==1)
180	{
181	for (j=1;j<=n;j++)
182	{
183	if (flaglist[j]==0){aux=aux+Exp[i][1][j];}
184	}
185	ss=aux; aux=0;
186	} // monomial
187	else
188	{
189	for (j=1;j<=n;j++)
190	{
191	if (flaglist[j]==0)
192	{
193	aux1=aux1+Exp[i][1][j];
194	aux2=aux2+Exp[i][2][j];
195	}
196	}
197	ss=aux1,aux2; aux1=0; aux2=0;
198	} // binomial
199	sums[i]=ss;}
200	return(sums);
201	}
202	example
203	{"EXAMPLE:"; echo = 2;
204	ring r = 0,(x(1),y(2),x(3),y(4),x(5..7),y(8)),dp;
205	list flag=identifyvar();
206	ideal J=x(1)^3x(3)-y(2)y(4)^2,x(5)y(2)-x(7)y(4)^2,x(6)^2(1-y(4)y(8)^5);
207	list L=data(J,3,8);
208	list EL=Edatalist(L[1],L[2],3,8,flag);
209	EL; // E-order of each term
210
211
212	ring r = 2,(x(1),y(2),x(3),y(4),x(5..7),y(8)),dp;
213	list flag=identifyvar();
214	ideal J=x(1)^3x(3)-y(2)y(4)^2,x(5)y(2)-x(7)y(4)^2,x(6)^2(1-y(4)y(8)^5);
215	list L=data(J,3,8);
216	list EL=Edatalist(L[1],L[2],3,8,flag);
217	EL; // E-order of each term IN CHAR 2, COMPUTATIONS NEED TO BE DONE IN CHAR 0
218
219
220	ring r = 0,(x(1..3)),dp;
221	list flag=identifyvar();
222	ideal J=x(1)^4*x(2)^2, x(1)^2-x(3)^3;
223	list L=data(J,2,3);
224	list EL=Edatalist(L[1],L[2],2,3,flag);
225	EL; // E-order of each term
226	}
227	///////////////////////////////////////////////////////////////////////////////////
228
229	proc EOrdlist(list Coef,list Exp,int k,int n,list flaglist)
230	"USAGE: EOrdlist(Coef,Exp,k,n,flaglist);
231	Coef,Exp,flaglist lists, k,n, integers
232	Exp is a list of lists of exponents, k=size(Exp)
233	COMPUTE: computes de E-order of an ideal given by a list (Coef,Exp) and extra information
234	RETURN: maximal E-order, and its position=number of generator and term
235	EXAMPLE: example EOrdlist; shows an example
236	"
237	{int i,can,canpost,lon;
238	number canmin;
239	list sums;
240
241	sums=Edatalist(Coef,Exp,k,n,flaglist);
242
243	canmin=sums[1][1]; // inicializating, works also with a monomial
244	for (i=1;i<=k; i++)
245	{
246	lon=size(sums[i]); // this is 2 for binomial and 1 for monomial generators
247	if (sums[i][1]<=canmin and Coef[i][1]!=0)
248	{
249	canmin=sums[i][1];
250	can=i; canpost=1;
251	}
252
253	// if the generator is a binomial we check the second term
254
255	if (lon==2)
256	{
257	if (sums[i][2]<canmin and Coef[i][2]!=0)
258	{
259	canmin=sums[i][2];
260	can=i; canpost=2;
261	}
262	}
263	}
264	return(canmin,can,canpost);
265	}
266	example
267	{"EXAMPLE:"; echo = 2;
268	ring r = 0,(x(1),y(2),x(3),y(4),x(5..7),y(8)),dp;
269	list flag=identifyvar();
270	ideal J=x(1)^3x(3)-y(2)y(4)^2,x(5)y(2)-x(7)y(4)^2,x(6)^2(1-y(4)y(8)^5),x(7)^4*y(8)^2;
271	list L=data(J,4,8);
272	list Eo=EOrdlist(L[1],L[2],4,8,flag);
273	Eo[1]; // E-order
274	Eo[2]; // generator giving the E-order
275	Eo[3]; // term giving the E-order
276	}
277
278	//////////////////////////////////////////////////////
279
280	proc maxEord(list Coef,list Exp,int k,int n,list flaglist)
281	"USAGE: maxEord(Coef,Exp,k,n,flaglist);
282	Coef,Exp,flaglist lists, k,n, integers
283	Exp is a list of lists of exponents, k=size(Exp)
284	RETURN: computes de maximal E-order of an ideal given by Coef,Exp
285	EXAMPLE: example maxEord; shows an example
286	"
287	{
288	int i,lon;
289	number canmin; // THE ASSIGNMENT IS NOT OK BECAUSE IT IS OF TYPE NUMBER
290	list sums;
291
292	sums=Edatalist(Coef,Exp,k,n,flaglist);
293
294	canmin=sums[1][1]; // inicializating, works also with a monomial
295	for (i=1;i<=k; i++)
296	{
297	lon=size(sums[i]); // this is 2 for binomial and 1 for monomial generators
298	if (sums[i][1]<=canmin and Coef[i][1]!=0)
299	{
300	canmin=sums[i][1];
301	}
302
303	// if the generator is a binomial we check the second term
304
305	if (lon==2)
306	{
307	if (sums[i][2]<canmin and Coef[i][2]!=0)
308	{
309	canmin=sums[i][2];
310	}
311	}
312	}
313	return(canmin,sums);
314	}
315	example
316	{"EXAMPLE:"; echo = 2;
317	ring r = 0,(x(1),y(2),x(3),y(4),x(5..7),y(8)),dp;
318	list flag=identifyvar();
319	ideal J=x(1)^3x(3)-y(2)y(4)^2x(3),x(5)y(2)-x(7)y(4)^2,x(6)^2(1-y(4)y(8)^5),x(7)^4y(8)^2;
320	list L=data(J,4,8);
321	list M=maxEord(L[1],L[2],4,8,flag);
322	M[1]; // E-order
323	}
324	//////////////////////////////////////////////////////
325
326	proc elimrep(list maxvar)
327	"USAGE: elimrep(L); L is a list
328	COMPUTE: Eliminate repeated terms from a list
329	RETURN: the same list without repeated terms
330	EXAMPLE: example elimrep; shows an example
331	"
332	{
333	int i,j;
334	list aux2;
335
336	aux2=maxvar;
337	for (i=1;i<=size(aux2); i++)
338	{
339	for (j=i+1;j<=size(aux2); j++)
340	{
341	if (aux2[i]==aux2[j] and i!=j){aux2=delete(aux2,j);}
342	}
343	}
344	maxvar=aux2;
345	return(maxvar);
346	}
347	example
348	{"EXAMPLE:"; echo = 2;
349	ring r = 0,(x(1..3)),dp;
350	list L=4,5,2,5,7,8,6,3,2;
351	elimrep(L);
352	}
353
354	proc cleanunit(list mon,int n,list flaglist)
355	"USAGE: cleanunit(mon,n,flaglist);
356	mon, flaglist lists, n integer
357	COMPUTE: We clean (or forget) the units in a monomial, given by "y" variables
358	RETURN: The list defining the monomial ideal already cleaned
359	EXAMPLE: example cleanunit; shows an example
360	"
361	{
362	int i;
363	for (i=1;i<=n;i++)
364	{
365	if (flaglist[i]==1){mon[i]=0;}
366	}
367	// coef[1]=coef[1]*y(i)^mon[i]; IS NOT ALLOWED because mon[i] can be a number
368	// therefore, the coefficients remain constant
369	return(mon);
370	}
371	example
372	{"EXAMPLE:"; echo = 2;
373	ring r = 0,(x(1),y(2),x(3),y(4)),dp;
374	list flag=identifyvar();
375	ideal J=x(1)^3y(2)x(3)^5*y(4)^8;
376	list L=data(J,1,4);
377	L[2][1][1]; // list of exponents of the monomial J
378	list M=cleanunit(L[2][1][1],4,flag);
379	M; // new list without units
380	}
381	//////////////////////////////////////////////////////
382	// Classification of the ideal E-Coeff_V(P):
383	// ccase=1, E-Coeff_V(P)=0
384	// 2,3 Bold regular case
385	// 4 P=1 monomial case (detected before)
386	// 0 Otherwise
387
388	proc ECoef(list Coef,list expP,int sP,int V,number auxc,int n,list flaglist)
389	"USAGE: ECoef(Coef,expP,sP,V,auxc,n,flaglist);
390	Coef, expP, flaglist lists, sP, V, n integers, auxc number
391	COMPUTE: The ideal E-Coeff_V(P), where V is a permissible hypersurface which belongs to the center
392	RETURN: list of exponents, list of coefficients and classification of the ideal E-Coeff_V(P)
393	EXAMPLE: example ECoef; shows an example
394	"
395	{
396	int i,j,k,l,numg,ccase,cont2,cont3,val;
397	number aa;
398	list Eco,newcoef,auxexp,newL,rs,rs2,aux,aux2,aux3,aux4,L;
399
400	auxexp=expP;
401
402	l=1;
403	for (i=1;i<=sP;i++)
404	{
405	rs[i]=size(Coef[i]);
406	if (rs[i]==2)
407	{ // binomials
408	if (auxexp[i][1][V]!=auxexp[i][2][V]) // no common factors for the variable in V
409
410	{
411	for (j=1;j<=2;j++)
412	{
413	if (auxexp[i][j][V]<auxc)
414	{
415	aa=auxc/(auxc-auxexp[i][j][V]);
416	auxexp[i][j][V]=0;
417	aux4[1]=multiplylist(auxexp[i][j],aa);
418	Eco[l]=aux4;
419	// newcoef[l]=Coef[i][j]^aa; IT IS NO ALLOWED!!!
420	newcoef[l]=Coef[i][j]; // we leave it constant
421	l=l+1;
422	}
423	}
424	}
425	else // common factors for the variable in V, of zero in both terms
426
427	{
428	if (auxexp[i][1][V]<auxc)
429	{
430	aa=auxc/(auxc-auxexp[i][1][V]);
431	auxexp[i][1][V]=0; auxexp[i][2][V]=0;
432
433	// this generator is a power of a binomial
434	// one possibility is Eco[l]=auxexp[i]; we leave it constant and add some extra number aa, or
435	// define a binomial again. The E-order coincides!!!
436
437	aux=multiplylist(auxexp[i][1],aa);
438	aux2=multiplylist(auxexp[i][2],aa);
439	aux3[1]=aux;
440	aux3[2]=aux2;
441	Eco[l]=aux3; newcoef[l]=Coef[i];
442	l=l+1;
443	}
444	}
445	}
446	else // monomials
447	{
448	if (auxexp[i][1][V]<auxc)
449	{
450	aa=auxc/(auxc-auxexp[i][1][V]);
451	auxexp[i][1][V]=0;
452	aux4=list();
453	aux4[1]=multiplylist(auxexp[i][1],aa);
454	Eco[l]=aux4;
455	newcoef[l]=Coef[i];
456	l=l+1;
457	}
458	}
459	}
460
461	// cleaning units from the monomial generators of Eco
462	// If there are hyperbolic equations in Eco, such that Eco=1, we detect it later, computing the E-order
463
464	numg=size(Eco);
465	for (k=1;k<=numg;k++)
466	{
467	if (size(newcoef[k])==1){Eco[k][1]=cleanunit(Eco[k][1],n,flaglist);}
468	}
469
470	// checking Eco
471
472	ccase=0;
473	cont2=0;
474	cont3=0;
475	val=0;
476
477	// CASE Eco=0: If Eco=empty list then as ideal Eco=0
478
479	if (numg==0){ccase=1;}
480	else
481	{
482	for (i=1;i<=numg;i++)
483	{
484	rs2[i]=size(newcoef[i]);
485	if (rs2[i]==1)
486	{
487	val=val+n; // monomials
488	for (l=1;l<=n; l++)
489	{
490	if (Eco[i][1][l]==0) {cont2=cont2+1;}
491	}
492	}
493	else
494	{
495	val=val+(2*n); // binomials
496	for (l=1;l<=n; l++)
497	{
498	if (Eco[i][1][l]==0) {cont2=cont2+1;}
499	if (Eco[i][2][l]==0) {cont2=cont2+1;}
500	}
501	}
502	}
503
504	// If cont2=val then all the entries of Eco are zero!! As ideal Eco=1
505
506	for (i=1;i<=sP;i++)
507	{
508	if (rs[i]==2)
509	{ // binomials
510	for (l=1;l<=n;l++)
511	{
512	if (expP[i][1][l]!=0) {cont3=cont3+1;}
513	if (expP[i][2][l]!=0) {cont3=cont3+1;}
514	}
515	}
516	else
517	{ // monomials
518	for (l=1;l<=n;l++)
519	{
520	if (expP[i][1][l]!=0)
521	{
522	cont3=cont3+1;
523	}
524	}
525	}
526	}
527
528	// If cont3=0 all the entries of expP are zero!! As ideal P=1 this is detected before
529	// If cont3=1 then P is bold regular
530
531
532	// CASE Eco=1
533
534	if (cont2==val and cont3==1){ccase=2;} // BOLD REGULAR CASE
535	if (cont2==val and cont3>1){ccase=3;} // CASE P=x^{\alpha},x^{\beta}, IN FACT, BOLD REGULAR
536	if (cont2==val and cont3==0){ccase=4;} // P=1, then I=1 monomial case
537
538	// Case BOLD REGULAR P=x^{\alpha}*(1-\mu y^{\delta})
539	// IT IS NON NECESSARY TO CHECK IT, Eco=empty list, already done!
540
541	L=maxEord(newcoef,Eco,numg,n,flaglist); // L[1] is the E-order of Eco
542	if (L[1]==0){ccase=2; print("E-order zero!");} // BOLD REGULAR CASE
543
544	// we leave it to check the computations
545
546	} // close else
547
548	return(Eco,newcoef,ccase);
549	}
550	example
551	{"EXAMPLE:"; echo = 2;
552	ring r = 0,(x(1),y(2),x(3),y(4),x(5..7)),dp;
553	list flag=identifyvar();
554	ideal P=x(1)^2x(3)^5-x(5)^7y(4),x(6)^3y(2)^5-x(7)^5,x(5)^3x(6)-y(4)^3*x(1)^5;
555	list L=data(P,3,7);
556	list L2=ECoef(L[1],L[2],3,1,3,7,flag);
557	L2[1]; // exponents of the E-Coefficient ideal respect to x(1)
558	L2[2]; // its coefficients
559	L2[3]; // classify the type of ideal obtained
560
561	ring r = 0,(x(1),y(2),x(3),y(4)),dp;
562	list flag=identifyvar();
563	ideal J=x(1)^3(1-2y(2)*y(4)^2); // Bold regular case
564	list L=data(J,1,4);
565	list L2=ECoef(L[1],L[2],1,1,3,4,flag);
566	L2;
567
568	ring r = 0,(x(1),y(2),x(3),y(4),x(5..7)),dp;
569	list flag=identifyvar();
570	ideal J=x(1)^3-x(3)^2y(4)^2,x(1)x(7)y(2)-x(6)^3x(5)y(4)^3,x(5)^3-x(5)^3y(2)^2;
571	list L=data(J,3,7);
572	list L2=ECoef(L[1],L[2],3,1,2,7,flag);
573	L2;
574
575	ring r = 3,(x(1),y(2),x(3),y(4),x(5..7)),dp;
576	list flag=identifyvar();
577	ideal J=x(1)^3-x(3)^2y(4)^2,x(1)x(7)y(2)-x(6)^3x(5)y(4)^3,x(5)^3-x(5)^3y(2)^2;
578	list L=data(J,3,7);
579	list L2=ECoef(L[1],L[2],3,1,2,7,flag);
580	L2; // THE COMPUTATIONS ARE NOT CORRECT IN CHARACTERISTIC p>0
581	// because numbers are treated as 0 in assignments
582
583	}
584	////////////////////////////////////////////////////////////////////////////
585	// The intvec a indicates the previous center
586	// Hhist = intvec of exceptional divisors of the parent chart
587
588	proc determinecenter(list Coef,list expJ,number c,int n,int Y,intvec a,list listmb,list flag,int control3,intvec Hhist)
589	"USAGE: determinecenter(Coef,expJ,c,n,Y,a,listmb,flag,control3,Hhist);
590	Coef, expJ, listmb, flag lists, c number, n, Y, control3 integers, a, Hhist intvec
591	COMPUTE: next center of blowing up and related information, see example
592	RETURN: several lists defining the center and related information
593	EXAMPLE: example determinecenter; shows an example
594	"
595	{
596	int i,j,rstep,l,mm,cont,cont1,cont2,cont3,a4,sI,sP,V,V2,ccase,b,Mindx,tip;
597	number auxc,a1,a2,ex,maxEo,aux;
598
599	list D,H,auxJ; // lists of D_n,D_n-1,...,D_1; H_n,H_n-1,...,H_1; J_n,J_n-1,...,J_1
600
601	list oldOlist,oldC,oldt,oldD,oldH,allH; // information of the previous step
602
603	list Olist,C,t,Dstar,center,expI,expP,newJ,maxset;
604
605	list maxvar,auxlist,aux3,auxD,auxolist,auxdiv,auxaux,L,rs,auxgamma,auxg2,aux1; // auxiliary lists
606	list auxinvlist,newcoef,EL,Ecoaux,Hplus,transH,Hsum,auxset,sumnewH; // auxiliary lists
607	list auxcoefI;
608
609	intvec oldinfobo7,infobo7;
610	int infaux,leh,leh2,leh3;
611
612	tip=listmb[1]; // It is not used in this procedure, it is used to compute the lcm of the denominators
613	oldOlist=listmb[2];
614	oldC=listmb[3];
615	oldt=listmb[4]; // t= resolution function
616	oldD=listmb[5];
617
618	oldH=listmb[6];
619	allH=listmb[7];
620
621	oldinfobo7=listmb[8]; // auxiliary intvec, it is used to define BO[7]
622
623	// inicializating lists
624	Olist=list();
625	C=list();
626	auxinvlist=list();
627
628	auxJ[1]=expJ;
629	rstep=n; // we are in dimension rstep
630	auxc=c;
631	cont=1;
632
633	if (Y==0) {D=iniD(n); H=iniD(n); infobo7=-1;} // first center, inicializate previous information
634
635	if (Y!=0 and rstep==n) // In dimension n, D'_n is always of this form
636	{ auxdiv=list0(n);
637	Dstar[1]=oldD[1];
638
639	b=size(a);
640	for (i=1;i<=n;i++) {for (j=1;j<=b;j++) {if (a[j]==i) {aux=aux+oldD[1][i];}}}
641	Dstar[1][Y]=aux;
642	aux=0;
643
644	auxdiv[Y]=oldOlist[1]-oldC[1];
645	D[1]=sumlist(Dstar[1],auxdiv);} // list defining D_n
646
647	// computing strict transforms of the exceptional divisors H
648
649	if (Y!=0){transH=iniD(n);
650	for (i=1;i<=size(oldH);i++){transH[i]=oldH[i]; transH[i][Y]=0;} // Note: size(oldH)<=n
651	allH[Y]=1;} // transform of \|H\|=H_nU...UH_1
652
653	// We put here size(oldH) instead of n because maybe we have not
654	// calculated all the dimensions in the previous step
655
656	// STARTING THE LOOP
657
658	while (rstep>=1)
659	{
660	if (Y!=0 and rstep!=n) // transformation law of D_i for i<n
661	{
662	if (cont!=0) // the resolution function did not drop in higher dimensions
663	{
664	if (oldt[n-rstep]==a1/a2 and c==oldC[1] and control3==0)
665	{auxD=list0(n);
666	auxD[Y]=oldOlist[n-rstep+1]-oldC[n-rstep+1];
667	Dstar[n-rstep+1]=oldD[n-rstep+1];
668
669	for (i=1;i<=n;i++) {for (j=1;j<=b;j++) {if (a[j]==i) {aux=aux+oldD[n-rstep+1][i];}}}
670	Dstar[n-rstep+1][Y]=aux;
671	aux=0;
672
673	D[n-rstep+1]=sumlist(Dstar[n-rstep+1],auxD);
674
675	}
676	else
677	{cont=0;
678	for (j=n-rstep+1;j<=n; j++){D[j]=list0(n);}
679	}
680	}
681	}
682
683	// Factorizing J=M*I
684
685	cont1=0;
686	for (i=1;i<=n;i++) {if (D[n-rstep+1][i]==0) {cont1=cont1+1;}} // if it fails write: listO(n)[i]
687
688	if (cont1==n) {expI=expJ;} // D[n-rstep+1]=0 (is a list of zeros)
689	else {
690	for (i=1;i<=size(expJ);i++)
691	{rs[i]=size(Coef[i]);
692	if (rs[i]==2){ aux1=list();
693	aux1[1]=reslist(expJ[i][1],D[n-rstep+1]);
694	aux1[2]=reslist(expJ[i][2],D[n-rstep+1]);
695	expI[i]=aux1;} // binomial
696	else {aux1=list();
697	aux1[1]=reslist(expJ[i][1],D[n-rstep+1]);
698	expI[i]=aux1;}} // monomial
699	}
700
701	// NOTE: coeficients of I = coeficients of J, because I and J differ in a monomial
702
703	// Detecting errors, negative exponents in expI
704
705	sI=size(expI);
706
707	for (i=1;i<=sI;i++)
708	{
709	rs[i]=size(Coef[i]);
710	if (rs[i]==2)
711	{
712	for (j=1;j<=2;j++)i
713	{
714	for (l=1;l<=n; l++)
715	{
716	if (expI[i][j][l]<0)
717	{
718	ERROR("the BINOMIAL ideal I has negative components");
719	}
720	}
721	}
722	}
723	else
724	{
725	for (l=1;l<=n; l++)
726	{
727	if (expI[i][1][l]<0)
728	{
729	"the MONOMIAL ideal I has negative components";
730	"M ideal";
731	print(D[n-rstep+1]); print(expI);
732	"dimension"; print(rstep);
733	}
734	}
735	}
736	}
737
738	// Compute the maximal E-order of I
739
740	L=maxEord(Coef,expI,sI,n,flag);
741	maxEo=L[1]; // E-order of I
742
743	// Inicializating information
744
745	auxolist=maxEo;
746	a1=maxEo;
747	a2=auxc;
748	Olist=Olist+auxolist; // list of new maximal E-orders o_n,o_{n-1},...o_1
749	aux3=auxc;
750	C=C+aux3; // list of new critical values c=c_{n+1},c_{n},...c_2
751
752	// It is necessary to check if the first coordinate of the invariant has dropped or not
753	// NOTE: By construction, the first coordinate is always 1 !!
754	// It has dropped is equivalent to: CURRENT C<c of the previous step
755
756	// Calculate new H, this is done for every dimension
757
758	if (Y!=0){a4=size(oldt);
759	if (n-rstep+1>a4){cont=0; oldt[n-rstep+1]=0; } // VERIFICAR!!!!
760
761	if (cont!=0 and oldt[n-rstep+1]==a1/a2 and c==oldC[1] and control3==0){H[n-rstep+1]=transH[n-rstep+1];
762
763	// we fill now the value for BO[7]
764	if (oldinfobo7[n-rstep+1]==-1){leh=size(Hhist);
765	infobo7[n-rstep+1]=Hhist[leh];} // suitable index !!!
766	else{ infaux=oldinfobo7[n-rstep+1];
767	infobo7[n-rstep+1]=infaux;} // the same as the previous step
768
769	}
770	else {
771	if (rstep<n) {sumnewH=list0(n);
772	for (i=1;i<n-rstep+1;i++){sumnewH=sumlist(sumnewH,H[i]);}
773	H[n-rstep+1]=reslist(allH,sumnewH);}
774	else {H[n-rstep+1]=allH;}
775
776	// we fill the value for BO[7] too, we complete it at the end if necessary
777	infobo7[n-rstep+1]=-1;
778	}
779	}
780
781	// It is necessary to detect the monomial case AFTER inicializate the information
782	// OTHERWISE WE WILL HAVE EMPTY COMPONENTS IN THE RESOLUTION FUNCTION
783
784	// If maxEo=0 but maxo!=0 MONOMIAL CASE (because E-Sing(J,c) still !=emptyset)
785	// If maxEo=0 and maxo=0 then I=1, (real) monomial case, the same case for us
786	// NOTE THAT IT DOESN'T MATTER IF THERE IS A p-TH POWER OF A HYPERBOLIC EQ, THE E-ORDER IS ZERO ANYWAY
787
788	if (maxEo==0){auxgamma=Gamma(D[n-rstep+1],auxc,n); // Gamma gives (maxlist,gamma,center)
789	auxg2=auxgamma[3];
790	center=center+auxg2;
791	center=elimrep(center);
792	auxinvlist=auxgamma[2]; break;}
793
794	// Calculate P // P=I+M^{o/(c-o)} with weight o
795
796	if (maxEo>=auxc) {expP=expI; Mindx=0;} // The coefficients also remain constant
797	else {ex=maxEo/(auxc-maxEo);
798	auxlist=list();
799	Mindx=1;
800	auxlist[1]=multiplylist(D[n-rstep+1],ex); // weighted monomial part: D[n-rstep+1]^ex;
801	expP=insert(expI,auxlist); // P=I+D[n-rstep+1]^ex;
802	auxcoefI=Coef;
803	Coef=insert(Coef,list(1));} // Adding the coefficient for M
804
805	// NOTE: IT IS NECESSARY TO ADD COEFFICIENT 1 TO THE MONOMIAL PART M
806	// E-ord(P_i)=E-ord(I_i) so to compute the E-order of P_i we can compute E-ord(I_i)
807
808	// Calculate variables of E-maximal contact, ALWAYS WITH RESPECT TO THE IDEAL I !!
809
810	sP=size(expP); // Can be different from size(expI)
811
812	if (Mindx==1){ maxvar=Emaxcont(auxcoefI,expI,sI,n,flag);}
813	else{ maxvar=Emaxcont(Coef,expP,sP,n,flag);}
814
815	auxc=maxvar[1]; // E-order of P, critical value for the next step, ALSO VALID auxc=maxEo;
816	if (auxc!=maxEo)
817	{
818	ERROR("the E-order is not well computed");
819	}
820
821	maxset=maxvar[2];
822	center=center + maxset;
823
824	// Cleaning the center: eliminating repeated variables
825
826	center=elimrep(center);
827
828	if (rstep==1) {break;} // Induction finished, is not necessary to compute the rest
829
830	// Calculate Hplus=set of non permissible hypersurfaces
831	// RESET Hplus if c has dropped or we have eliminated hyperbolic generators
832
833	// ES NECESARIO PONER CONDICION DE SI INVARIANTE BAJA O NO??? SI BAJA HPLUS NO SE USA...
834
835	if (Y==0 or c<oldC[1] or control3==1) {Hplus=list0(n);}
836	else {Hsum=list0(n);
837	Hplus=allH;
838	for (i=1;i<=n-rstep+1;i++){Hsum=sumlist(Hsum,H[i]);}
839	Hplus=reslist(Hplus,Hsum); // CHEQUEAR QUE NO SALEN -1'S
840	}
841
842	// Taking into account variables of maxset outside of Hplus (so inside Hminus)
843
844	if (Y==0 or c<oldC[1] or control3==1){V=maxset[1]; // Hplus=0 so any variable is permissible
845	maxset=delete(maxset,1);} // eliminating this variable V from maxset
846	else{
847	// If the invariant remains constant V comes from the previous step
848
849	if (cont!=0 and oldt[n-rstep+1]==a1/a2 and c==oldC[1]){
850	if (Mindx==1){
851	//----------------------------USING HPLUS----------------------------------------
852	// REMIND THAT IN THIS CASE maxset=HYPERSURFACES OF E-MAXIMAL CONTACT FOR I, INSTEAD OF P
853
854	cont2=1;
855	cont3=1;
856	auxset=maxset;
857	while (cont2!=0){mm=auxset[1];
858	if (Hplus[mm]!=0) {auxset=delete(auxset,1); cont3=cont3+1;}
859	// eliminating non permissible variables from maxset
860	else {cont2=0;}}
861	V=maxset[cont3]; // first permissible variable
862	maxset=delete(maxset,cont3);
863	V2=a[n-rstep+1]; // can be different from the variable coming from the previous step
864	}
865
866	//-------------------------------------------------------------------------------
867	else{ V=a[n-rstep+1];}
868	}
869	else {V=maxset[1]; // Hplus=0 so any variable is permissible
870	maxset=delete(maxset,1);
871	}
872
873	}
874
875	// Calculate Eco=E-Coeff_V(P) where V is a permissible hypersurface which belongs to the center
876	// Eco can have rational exponents
877
878	Ecoaux=ECoef(Coef,expP,sP,V,auxc,n,flag);
879
880	// SPECIAL CASES: BOLD REGULAR CASE
881	//--------------------------------------------------------------------
882
883	if (Ecoaux[3]==1)
884	{ // Eco=EMPTY LIST, Eco=0 AS IDEAL
885	aux1[1]=list0(n);
886	newJ[1]=aux1; // monomial with zero entries, newJ=1 as ideal
887	newcoef[1]=list(1); // the new coefficient is only 1
888	auxaux=list();
889	auxaux[1]=newJ;
890	auxJ=auxJ+auxaux; // auxJ list of ideals J_i
891	auxinvlist=1;
892	break;
893	}
894	//-----------------------------------------------------------
895	// THIS CASE IS NOT GOING TO APPEAR, BUT WE LEAVE IT TO CHECK COMPUTATIONS
896	if (Ecoaux[3]==2 or Ecoaux[3]==3)
897	{ // Eco=0 LIST, Eco=1 AS IDEAL
898	aux1[1]=list0(n);
899	newJ[1]=aux1;
900	newcoef[1]=list(1); print("Strange case happens");
901	auxaux=list();
902	auxaux[1]=newJ;
903	auxJ=auxJ + auxaux; // auxJ list of ideals J_i
904	auxinvlist=1;
905	break;
906	}
907	//-----------------------------------------------------------
908	// THIS CASE IS NOT GOING TO APPEAR, BUT WE LEAVE IT TO CHECK COMPUTATIONS
909
910	// P=1 THIS CANNOT HAPPEN SINCE P=1 IFF I=1 (or I is equivalent to 1)
911	// and this is the monomial case, already checked
912
913	if (Ecoaux[3]==4){ERROR("ERROR in ECoef");}
914	//-----------------------------------------------------------
915	// If we are here Ecoaux[3]=0, then continue
916	// Filling the list of "ideals", auxJ=J_n,J_{n-1},...
917	newJ=Ecoaux[1];
918	newcoef=Ecoaux[2];
919
920	auxJ=insert(auxJ,newJ,n-rstep+1); // newJ is inserted after n-rstep+1 position, so in position n-rstep+2
921
922	// New input for the loop, if we are here newJ is different from 0
923
924	expJ=newJ;
925	Coef=newcoef;
926
927	newJ=list();
928	expI=list();
929	expP=list();
930	rstep=rstep-1; // print(size(auxJ));
931	}
932	// EXIT LOOP "while"
933	// we do NOT construct the center as an ideal because WE USE LISTS
934	t=dividelist(Olist,C); // resolution function t
935	// Complete the intvec infobo7 if necessary
936	if (control3==1){infobo7=-1;} // We reset the value after clean hyperbolic equations
937	leh2=size(Olist);
938	leh3=size(infobo7);
939	if (leh3<leh2){for (j=leh3+1;j<=leh2; j++){infobo7[j]=-1;}}
940	// Auxiliary list to complete the resolution function in special cases
941	if (size(auxinvlist)==0) {auxinvlist[1]=0;}
942	return(center,auxJ,Olist,C,t,D,H,allH,auxinvlist,infobo7);
943	}
944	example
945	{"EXAMPLE:"; echo = 2;
946	ring r = 0,(x(1..4)),dp;
947	list flag=identifyvar();
948	ideal J=x(1)^2-x(2)^2x(3)^5, x(1)x(3)^3+x(4)^6;
949	list Lmb=1,list0(4),list0(4),list0(4),iniD(4),iniD(4),list0(4),-1;
950	list L=data(J,2,4);
951	list LL=determinecenter(L[1],L[2],2,4,0,0,Lmb,flag,0,-1); // Compute the first center
952	LL[1]; // index of variables in the center
953	LL[2]; // exponents of ideals J_4,J_3,J_2,J_1
954	LL[3]; // list of orders of J_4,J_3,J_2,J_1
955	LL[4]; // list of critical values
956	LL[5]; // components of the resolution function t
957	LL[6]; // list of D_4,D_3,D_2,D_1
958	LL[7]; // list of H_4,H_3,H_2,H_1 (exceptional divisors)
959	LL[8]; // list of all exceptional divisors acumulated
960	LL[9]; // auxiliary invariant
961	LL[10]; // intvec pointing out the last step where the function t has dropped
962
963	ring r= 0,(x(1..4)),dp;
964	list flag=identifyvar();
965	ideal J=x(1)^3-x(2)^2x(3)^5, x(1)x(3)^3+x(4)^5;
966	list Lmb=2,list0(4),list0(4),list0(4),iniD(4),iniD(4),list0(4),-1;
967	list L2=data(J,2,4);
968	list L3=determinecenter(L2[1],L2[2],2,4,0,0,Lmb,flag,0,-1); // Example with rational exponents in E-Coeff
969	L3[1]; // index of variables in the center
970	L3[2]; // exponents of ideals J_4,J_3,J_2,J_1
971	L3[3]; // list of orders of J_4,J_3,J_2,J_1
972	L3[4]; // list of critical values
973	L3[5]; // components of the resolution function
974	}
975	////////////////////////////////////////////////////////
976	// idchart= identity number of the current chart
977	// infochart=chart[idchart] information related to the chart to blow up
978	// infochart= int parent,int Y,intvec a,list expJ,list Coef, list flag, // NEEDED FOR THE RESOLUTION
979	// intvec Hhist, list blwhist, module path, list hipercoef, list hiperexp // NEEDED FOR THE OUTPUT
980
981	// NOTE: IT IS NOT NECESSARY TAKE INTO ACCOUNT "y" VARIABLES BECAUSE THE CENTER IS ALREADY GIVEN
982
983	proc Blowupcenter(list center,int idchart,int nchart,list infochart,number c,int n,int currentstep)
984	"USAGE: Blowupcenter(center,id,nchart,infochart,c,n,cstep);
985	center, infochart lists, id, nchart, n, cstep integers, c number
986	COMPUTE: The blowing up at the chart IDCHART along the given center
987	RETURN: new affine charts and related information, see example
988	EXAMPLE: example Blowupcenter; shows an example
989	"
990	{int num,i,j,k,l,parent,Y,lon,m,m2;
991	intvec a,Hhist,auxHhist;
992	number auxsum, auxsum2;
993	list sons,aux,expJ,blexpJ,blD;
994	list auxstep,Coef;
995	list auxchart,auxchart1,info,flaglist;
996	list auxblwhist,blwhist,hipercoef,hiperexp;
997	module auxpath,auxp2;
998
999	parent=idchart;
1000	num=size(center);
1001
1002	// Transform to intvec the list of variables defining the center
1003	a=center[1];
1004	for (i=2;i<=num;i++){a=a,center[i];}
1005
1006	expJ=infochart[4];
1007	Coef=infochart[5];
1008	flaglist=infochart[6];
1009	Hhist=infochart[7];
1010	blwhist=infochart[8];
1011	auxpath=infochart[9];
1012	hipercoef=infochart[10];
1013	hiperexp=infochart[11];
1014
1015	l=size(expJ);
1016
1017	// input for the loop
1018	blexpJ=expJ;
1019
1020	// making the blowing up in the i-th chart
1021	for (i=1;i<=num;i++)
1022	{
1023	// we assign the current number of charts +1 to the i-th chart
1024	idchart=nchart+1;
1025	nchart=nchart+1;
1026	aux=idchart;
1027	sons=sons+aux;
1028
1029	auxstep[i]=currentstep+1;
1030
1031	Y=center[i];
1032
1033	// The blowing up
1034
1035	for (j=1;j<=l;j++){lon=size(Coef[j]);
1036	if (lon==1){for (m=1;m<=n;m++){for (m2=1;m2<=num;m2++){
1037	if (m==center[m2]){auxsum=auxsum+ expJ[j][1][m];}}}
1038	blexpJ[j][1][Y]=auxsum-c;
1039	auxsum=0;} // monomial
1040	else {for (m=1;m<=n;m++){for (m2=1;m2<=num;m2++){
1041	if (m==center[m2]){auxsum=auxsum+expJ[j][1][m];
1042	auxsum2=auxsum2+expJ[j][2][m];}}}
1043	blexpJ[j][1][Y]=auxsum-c;
1044	blexpJ[j][2][Y]=auxsum2-c;
1045	auxsum=0; auxsum2=0;} // binomial
1046	}
1047
1048
1049	auxHhist=Hhist,Y; // history of the exceptional divisors in this chart
1050	auxblwhist=tradblwup(blwhist,n,Y,a,num); // history of the blow ups in this chart
1051
1052	auxp2=auxpath,[parent,i];
1053
1054	auxchart1=parent,Y,a,blexpJ,Coef,flaglist,auxHhist,auxblwhist,auxp2,hipercoef,hiperexp;
1055
1056	// Coef, flaglist are not modified after the blowing-up, the hyperbolic information is the same as in the parent chart
1057
1058	auxchart[i]=auxchart1;
1059
1060	// Inicializating the exponents of J for the next chart
1061
1062	blexpJ=expJ;
1063	}
1064	// end of the loop
1065
1066	// we add its sons to the current chart
1067	infochart=infochart+sons;
1068	info[1]=infochart;
1069
1070	return(info,auxchart,nchart,auxstep,num);
1071	}
1072	example
1073	{"EXAMPLE:"; echo = 2;
1074	ring r = 0,(x(1),y(2),x(3),y(4),x(5..7)),dp;
1075	list flag=identifyvar();
1076	ideal J=x(1)^3-x(3)^2y(4)^2,x(1)x(7)y(2)-x(6)^3x(5)y(4)^3,x(5)^3-x(5)^3y(2)^2;
1077	list Lmb=2,list0(7),list0(7),list0(7),iniD(7),iniD(7),list0(7),-1;
1078	list L=data(J,3,7);
1079	list L2=determinecenter(L[1],L[2],2,7,0,0,Lmb,flag,0,-1); // Computing the center
1080	module auxpath=[0,-1];
1081	list infochart=0,0,0,L[2],L[1],flag,0,list0(7),auxpath,list(),list();
1082	list L3=Blowupcenter(L2[1],1,1,infochart,2,7,0);
1083	L3[1]; // current chart (parent,Y,center,expJ,Coef,flag,Hhist,blwhist,path,hipercoef,hiperexp) with sons: [12],...,[16]
1084	L3[2][1]; // information of its first son, write L3[2][2],...,L3[2][5] to see the other sons
1085	L3[3]; // current number of charts
1086	L3[4]; // step/level associated to each son
1087	L3[5]; // number of variables in the center
1088	}
1089	//////////////////////////////////////////////////////////////
1090
1091	proc tradblwup(list blwhist,int n,int Y,intvec a,int num)
1092	"Internal procedure - no help and no example available
1093	"
1094	{
1095	int i,j,blwnew;
1096	intvec aux,aux2;
1097
1098	for (j=1;j<=n;j++){
1099	for (i=1;i<=num;i++){
1100	if (j==a[i] and a[i]!=Y){blwnew=Y; break;}
1101	else {blwnew=0;}
1102	}
1103	aux=blwhist[j];
1104	aux2=aux,blwnew;
1105	blwhist[j]=aux2;
1106	}
1107	return(blwhist);
1108	}
1109	//////////////////////////////////////////////////////////////
1110	// It is called only when Eord(J)=0, and J!=1 it is checked inside
1111	// SO IT IS CALLED AFTER: maxEord(Coef,expJ,sJ,n,flaglist); --> gives (max E-order,sums)
1112
1113	proc Nonhyp(list Coef,list expJ,int sJ,int n,list flaglist,list sums)
1114	"USAGE: Nonhyp(Coef,expJ,sJ,n,flaglist,sums);
1115	Coef, expJ, flaglist, sums lists, sJ, n integers
1116	COMPUTE: The "ideal" generated by the non hyperbolic generators of J
1117	RETURN: lists with the following information
1118	newcoef,newJ: coefficients and exponents of the non hyperbolic generators
1119	totalhyp,totalgen: coefficients and exponents of the hyperbolic generators
1120	flaglist: new list saying status of variables
1121	NOTE: the basering r is supposed to be a polynomial ring K[x,y],
1122	in fact, we work in a localization of K[x,y], of type K[x,y]_y with y invertible variables.
1123	EXAMPLE: example Nonhyp; shows an example
1124	"
1125	{
1126	int i,j,k,h,lon,lon2,cont;
1127	number eordcontrol;
1128	list genhyp,listgen,listid,posnumJ,newJ,newcoef,hypcoef,hyp,aux1,aux2,aux3,aux,midlist;
1129	list totalhyp,totalgen;
1130
1131	eordcontrol=0;
1132
1133	while (eordcontrol==0 and sJ!=0)
1134	{
1135
1136	// Give a positional number/flag to each generator of expJ
1137
1138	for (i=1;i<=sJ; i++){listgen=expJ[i]; listid=i; posnumJ[i]=listgen+listid; }
1139
1140	// Select the non hyperbolic and hyperbolic generators
1141
1142	for (j=1;j<=sJ; j++){lon=size(Coef[j]);
1143	if (lon==1){
1144
1145	// IS NOT NECESSARY TO CHECK IF THERE EXIST A MONOMIAL WITH ONLY UNITS, ALREADY DONE!!
1146
1147	aux1=aux1+posnumJ[j];
1148	aux3=list();
1149	aux3[1]=expJ[j];
1150	newJ=newJ+aux3;
1151	aux3[1]=Coef[j];
1152	newcoef=newcoef+aux3;
1153	}
1154
1155	else{ // CHECKING BINOMIALS, ONE TERM WITH E-ORDER ZERO GIVES HYPERBOLIC EQ
1156
1157	if (sums[j][1]==0 or sums[j][2]==0){aux2=aux2+posnumJ[j];
1158	aux3=list();
1159	aux3[1]=expJ[j];
1160	genhyp=genhyp+aux3;
1161	aux3[1]=Coef[j];
1162	hypcoef=hypcoef+aux3;
1163	if (sums[j][1]==0){aux3[1]=expJ[j][2]; hyp=hyp+aux3;}
1164	if (sums[j][2]==0){aux3[1]=expJ[j][1]; hyp=hyp+aux3;}
1165	}
1166	else {aux1=aux1+posnumJ[j];
1167	aux3=list();
1168	aux3[1]=expJ[j];
1169	newJ=newJ+aux3;
1170	aux3[1]=Coef[j];
1171	newcoef=newcoef+aux3;}
1172
1173	}
1174	}
1175
1176	// NOTE: aux1 and aux2 are no needed right now!
1177
1178	// Identify new y variables, that is, x variables in the monomials contained in hyp
1179
1180	h=size(hyp);
1181
1182	for (k=1;k<=h; k++){ for(i=1;i<=n; i++){ if (hyp[k][i]!=0 and flaglist[i]==0) {flaglist[i]=1;}}}
1183
1184	// To replace x by y IT IS NECESSARY TO CHANGE THE BASERING!!! We change only the list flaglist
1185
1186	// CHECK IF THE IDEAL IS ALREADY GENERATED BY MONOMIALS, in this case
1187	// WE HAVE FINISHED THE E-RESOLUTION PART, J GENERATED BY MONOMIALS AND HYPERBOLIC EQS
1188
1189	cont=0;
1190	lon2=size(newJ);
1191	for (j=1;j<=lon2; j++){if (size(newJ[j])==1){cont=cont+1;}}
1192
1193	if (cont==lon2){newcoef=list();
1194	newJ=list();
1195	totalgen=totalgen+genhyp;
1196	totalhyp=totalhyp+hypcoef;
1197	break;}
1198
1199	// CHECK IF THERE ARE MORE HYPERBOLIC EQUATIONS AFTER UPDATE THE FLAG LIST
1200	// CHECK THE MAXIMAL E-ORDER AGAIN
1201
1202	if (lon2==0){ // we are in the previous case, newJ=empty list, save values and exit
1203
1204	totalgen=totalgen+genhyp;
1205	totalhyp=totalhyp+hypcoef;
1206	break;
1207	}
1208
1209	midlist=maxEord(newcoef,newJ,lon2,n,flaglist);
1210
1211	eordcontrol=midlist[1];
1212
1213	if (eordcontrol==0){ // new input for the loop
1214	Coef=newcoef;
1215	expJ=newJ;
1216	sJ=lon2;
1217	sums=midlist[2]; // flaglist is already updated
1218
1219	totalgen=totalgen+genhyp;
1220	totalhyp=totalhyp+hypcoef;
1221
1222	hypcoef=list();
1223	genhyp=list();
1224
1225	newJ=list();
1226	newcoef=list();
1227	}
1228	else{ // If the process is already finished we save the values and exit
1229
1230	totalgen=totalgen+genhyp;
1231	totalhyp=totalhyp+hypcoef;
1232	}
1233
1234	} // closing while
1235
1236	return(newcoef,newJ,totalhyp,totalgen,flaglist);
1237	}
1238	example
1239	{"EXAMPLE:"; echo = 2;
1240	ring r = 0,(x(1),y(2),x(3),y(4),x(5..7)),dp;
1241	list flag=identifyvar(); // List giving flag=1 to invertible variables: y(2),y(4)
1242	ideal J=x(1)^3-x(3)^2y(4)^2,x(1)x(7)y(2)-x(6)^3x(5)y(4)^3,1-x(5)^2y(2)^2;
1243	list L=data(J,3,7);
1244	list L2=maxEord(L[1],L[2],3,7,flag);
1245	L2[1]; // Maximum E-order
1246	list New=Nonhyp(L[1],L[2],3,7,flag,L2[2]);
1247	New[1]; // Coefficients of the non hyperbolic part
1248	New[2]; // Exponents of the non hyperbolic part
1249	New[3]; // Coefficients of the hyperbolic part
1250	New[4]; // New hyperbolic equations
1251	New[5]; // New list giving flag=1 to invertible variables: y(2),y(4),y(5)
1252
1253	ring r = 0,(x(1..4)),dp;
1254	list flag=identifyvar();
1255	ideal J=1-x(1)^5x(2)^2x(3)^5, x(1)^2x(3)^3+x(1)^4x(4)^6;
1256	list L=data(J,2,4);
1257	list L2=maxEord(L[1],L[2],2,4,flag);
1258	L2[1]; // Maximum E-order
1259	list New=Nonhyp(L[1],L[2],2,4,flag,L2[2]);
1260	New;
1261
1262	}
1263	//////////////////////////////////////////////////////////////
1264
1265	proc calculateI(list Coef,list expJ,number c,int n,int Y,intvec a,number oldordI,list oldD)
1266	"USAGE: calculateI(Coef,expJ,c,n,Y,a,b,D);
1267	Coef, expJ, D lists, c, b numbers, n,Y integers, a intvec
1268	RETURN: ideal I, non monomial part of J
1269	EXAMPLE: example calculateI; shows an example
1270	"
1271	{
1272	int i,cont1,b,j;
1273	number EordI,aux;
1274	list D,L,expI;
1275	list auxdiv,Dstar,aux1,rs;
1276
1277	// WE NEED THE MONOMIAL PART, BUT ONLY IN DIMENSION n
1278
1279	auxdiv=list0(n);
1280	auxdiv[Y]=oldordI-c;
1281	Dstar[1]=oldD[1];
1282
1283	b=size(a);
1284	for (i=1;i<=n;i++) {for (j=1;j<=b;j++) {if (a[j]==i) {aux=aux+oldD[1][i];}}}
1285	Dstar[1][Y]=aux;
1286	aux=0;
1287
1288	D[1]=sumlist(Dstar[1],auxdiv);
1289
1290	cont1=0;
1291	for (i=1;i<=n;i++) {if (D[1][i]==0) {cont1=cont1+1;}} // if it fails write listO(n)[i]
1292
1293	if (cont1==n) {expI=expJ;}
1294	else {
1295	for (i=1;i<=size(expJ);i++)
1296	{rs[i]=size(Coef[i]);
1297	if (rs[i]==2){ aux1=list();
1298	aux1[1]=reslist(expJ[i][1],D[1]);
1299	aux1[2]=reslist(expJ[i][2],D[1]);
1300	expI[i]=aux1;} // binomial
1301	else {aux1=list();
1302	aux1[1]=reslist(expJ[i][1],D[1]);
1303	expI[i]=aux1;}} // monomial
1304	}
1305
1306	return(expI);
1307	}
1308	example
1309	{"EXAMPLE:"; echo = 2;
1310	ring r = 0,(x(1..3)),dp;
1311	list flag=identifyvar();
1312	ideal J=x(1)^4*x(2)^2, x(3)^3;
1313	list Lmb=1,list0(3),list0(3),list0(3),iniD(3),iniD(3),list0(3),-1;
1314	list L=data(J,2,3);
1315	list LL=determinecenter(L[1],L[2],3,3,0,0,Lmb,flag,0,-1); // Calculate the center
1316	module auxpath=[0,-1];
1317	list infochart=0,0,0,L[2],L[1],flag,0,list0(3),auxpath,list(),list();
1318	list L3=Blowupcenter(LL[1],1,1,infochart,3,3,0); // blowing-up and looking to the x(3) chart
1319	calculateI(L3[2][1][5],L3[2][1][4],3,3,3,L3[2][1][3],3,iniD(3)); // (I_3)
1320	// looking to the x(1) chart
1321	calculateI(L3[2][2][5],L3[2][2][4],3,3,1,L3[2][2][3],3,iniD(3)); // (I_3)
1322	}
1323	//////////////////////////////////////////////////////////////////////////////////////
1324	// //
1325	// E-RESOLUTION: Eresol(J) subroutine computing the E-resolution of J, char 0 //
1326	// //
1327	//////////////////////////////////////////////////////////////////////////////////////
1328
1329	proc Eresol(ideal J)
1330	"USAGE: Eresol(J); J ideal
1331	RETURN: The E-resolution of singularities of J in terms of the affine charts, see example
1332	EXAMPLE: example Eresol; shows an example
1333	"
1334	{int i,n,k,idchart,nchart,parent,Y,oldid,tnum,s,cont,control,control2,control3,cont2,val,rs2,l,cont3,tip;
1335	intvec a,Hhist;
1336	number c,EordJ,EordI,oldordI;
1337	list L,LL,oldD,t,auxL,finalchart,chart,auxchart,newL,auxp,auxfchart,L2;
1338	list Coef,expJ,expI,sons,oldOlist,oldC,oldt,oldH,allH,auxordJ,auxordI,auxmb,mobile,invariant;
1339	list step,nsons,auxinv,extraL,totalinv,auxsum;
1340	string empstring;
1341	module auxpath;
1342
1343	// ADDED LATER
1344
1345	list flag,newflag,blwhist,hipercoef,hiperexp,hipercoefson,hiperexpson;
1346	intvec infobo7;
1347
1348	export finalchart;
1349	// export nsons;
1350	// export tnum;
1351	// export nchart;
1352	// export step;
1353	export invariant;
1354	export auxinv;
1355	export mobile;
1356
1357	n=nvars(basering);
1358	flag=identifyvar();
1359
1360	k=size(J);
1361	// Checking input data
1362	if (inidata(J,k)==0){return("This library only works for binomial ideals.");}
1363
1364	idchart=1;
1365	nchart=1;
1366	parent=0;
1367	step=0;
1368	control=0;
1369	control2=0;
1370	control3=0;
1371
1372	// Translate the input ideal to a list
1373	auxL=data(J,k,n); // data gives (Coef,Exp)
1374
1375	// THEREAFTER WE WORK ALL THE TIME WITH LISTS
1376
1377	L=maxEord(auxL[1],auxL[2],k,n,flag); // gives (max E-ord, sums)
1378	EordJ=L[1];
1379
1380	// before the first blow up I=J
1381	EordI=EordJ;
1382
1383	// main loop AT EACH CHART WE MUST INICIALIZATE ALL THE VALUES AND
1384	// CONSTRUCT THE FIRST CHART chart[1] BEFORE THE LOOP
1385
1386	// at the first step, before the blow up, there are no exceptional divisors, Y=0
1387	Y=0;
1388	expJ=auxL[2];
1389	Coef=auxL[1];
1390	Hhist=0;
1391	blwhist=list0(n);
1392	auxpath=[0,-1];
1393	hipercoef=list(); // this is for the first chart
1394	hiperexp=list();
1395	auxp=parent,Y,a,expJ,Coef,flag,Hhist,blwhist,auxpath,hipercoef,hiperexp;
1396	chart[1]=auxp; // information of the first chart
1397
1398	tip=1;
1399	oldOlist=list0(n);
1400	oldC=list0(n);
1401	oldC[1]=EordJ; // non necessary here
1402	c=EordJ; // the value c is given by the previous step
1403	oldt=list0(n);
1404	oldD=iniD(n);
1405	oldH=iniD(n);
1406	allH=list0(n);
1407
1408	for (i=1;i<=n;i++){infobo7[i]=-1;}
1409
1410	auxmb=tip,oldOlist,oldC,oldt,oldD,oldH,allH,infobo7;
1411	mobile[1]=auxmb; // mobile corresponding to the first chart
1412	auxinv[1]=list(0);
1413
1414	// NOTE: oldC[1] is the value c to classify the chart in one of the next cases
1415
1416	// HERE BEGIN THE LOOP
1417
1418	while (idchart<=nchart) // WE PROCEED WHILE THERE EXIST UNSOLVED CHARTS
1419	{
1420	if (idchart!=1) // WE ARE NOT IN THE FIRST CHART, INICIALIZATE ALL THE VALUES
1421	{
1422
1423	parent=chart[idchart][1];
1424	Y=chart[idchart][2];
1425	a=chart[idchart][3];
1426	expJ=chart[idchart][4];
1427	Coef=chart[idchart][5];
1428	flag=chart[idchart][6];
1429	Hhist=chart[idchart][7]; // it is not necessary for the computations
1430	blwhist=chart[idchart][8];
1431	auxpath=chart[idchart][9];
1432	hipercoef=chart[idchart][10];
1433	hiperexp=chart[idchart][11];
1434
1435	k=size(Coef); // IT IS NECESSARY TO COMPUTE IT BECAUSE IT DECREASES IF THERE ARE HYPERBOLIC EQS
1436
1437	auxordJ=maxEord(Coef,expJ,k,n,flag);
1438	EordJ=auxordJ[1];
1439
1440	if (control==0){c=mobile[parent+1][3][1];} // we keep c from the last step
1441	else {c=EordJ; control=0; } // we reset the value of c
1442
1443	if (control2==1){c=EordJ; control2=0; control3=1;} // we reset the value of c
1444
1445	// NOTE: oldC[1] is the value c to classify the chart in one of the next cases
1446
1447	}
1448
1449	// The E-order must be computed here
1450
1451	oldid=idchart;
1452
1453	if (EordJ<0) {print("ERROR in J in chart"); print(idchart); break;}
1454
1455
1456	//-------------------------------------------------------------
1457	// CASE J=1, if we reset c, can happen Eord=c=0
1458
1459	// or if there are hyperbolic equations at the beginning!!!
1460
1461	// if (EordJ==0){auxfchart[1]=chart[idchart]; // WE HAVE FINISHED
1462	// finalchart=finalchart+auxfchart;
1463	// empstring="#"; print("reset c and Eord=c=0"); print(idchart);
1464	// invariant[idchart]=empstring;
1465	// auxinv[idchart]=list(0);
1466	// nsons[idchart]=0;
1467	// idchart=idchart+1;}
1468
1469
1470	//----------------------------------------------------------------------
1471	if (EordJ>=c and EordJ!=0) // subroutine: E-RESOLUTION OF PAIRS
1472	{
1473	if (parent>0)
1474	{ LL=determinecenter(Coef,expJ,c,n,Y,a,mobile[parent+1],flag,control3,chart[parent][7]); }
1475	else { LL=determinecenter(Coef,expJ,c,n,Y,a,mobile[parent+1],flag,control3,Hhist); }
1476
1477	// determinecenter gives (center,auxJ,Olist,C,t,D,H,allH,auxinvlist,infobo7)
1478
1479	// save current values, before the blow up
1480	oldOlist=LL[3];
1481	tip=computemcm(oldOlist);
1482	oldC=LL[4];
1483	oldt=LL[5];
1484	oldD=LL[6];
1485	oldH=LL[7];
1486	allH=LL[8];
1487	auxinv[idchart]=LL[9];
1488	infobo7=LL[10];
1489
1490	auxmb=tip,oldOlist,oldC,oldt,oldD,oldH,allH,infobo7;
1491	mobile[idchart+1]=auxmb;
1492	invariant[idchart]=oldt;
1493
1494	newL=Blowupcenter(LL[1],idchart,nchart,chart[idchart],c,n,step[idchart]);
1495
1496	// Blowupcenter gives (info,auxchart,nchart,auxstep,num)
1497
1498	// IMPORTANT: ADD THE NEW CHARTS AFTER EACH BLOW UP, IN ORDER TO KEEP THEM CORRECTLY
1499
1500	step=step+newL[4];
1501	nsons[idchart]=newL[5];
1502
1503	chart=chart+newL[2];
1504	finalchart=finalchart+newL[1];
1505
1506	// new input for the loop
1507
1508	idchart=idchart+1;
1509	nchart=newL[3];
1510
1511	control3=0;
1512
1513	} // END OF CASE EordJ>=c
1514	//---------------------------------------------------------------------
1515
1516	else{
1517
1518	// compute EordI=max E-order(I)
1519
1520	expI=calculateI(Coef,expJ,c,n,Y,a,mobile[parent+1][2][1],mobile[parent+1][5]);
1521	k=size(expJ); // probably non necessary
1522	auxordI=maxEord(Coef,expI,k,n,flag);
1523	EordI=auxordI[1];
1524	auxsum=auxordI[2];
1525
1526	// CASE EordI>0 DROP c AND CONTINUE
1527
1528	if (EordI>0){idchart=idchart; // keep the chart and back to the main loop while, dropping the value of c
1529	control=1;}
1530	else{ // EordI=0, so check if I=1 or not
1531
1532	cont2=0; // If cont2=val then all the entries of expI are zero!!
1533	val=0;
1534
1535	for (i=1;i<=k;i++) {rs2=size(Coef[i]);
1536	if (rs2==1){if (auxsum[i][1]==0){cont2=val; break;} // THERE EXIST A MONOMIAL WITH ONLY UNITS
1537
1538	val=val+n; // monomials
1539	for (l=1;l<=n; l++) {if (expI[i][1][l]==0) {cont2=cont2+1;}}
1540	}
1541	else{val=val+(2*n); // binomials
1542	for (l=1;l<=n; l++) {if (expI[i][1][l]==0) {cont2=cont2+1;}
1543	if (expI[i][2][l]==0) {cont2=cont2+1;}}
1544	}
1545	}
1546
1547
1548	// CASE EordI==0 AND I=1 THIS CHART IS DONE, FINISH
1549
1550	// NOTE: THIS CASE IS NOT MONOMIAL BECAUSE E-Sing(J,c) is empty
1551
1552	if (cont2==val){auxfchart[1]=chart[idchart];
1553	finalchart=finalchart+auxfchart;
1554	empstring="#";
1555	invariant[idchart]=empstring;
1556	auxinv[idchart]=list(0);
1557	nsons[idchart]=0;
1558
1559	// information for the mobile
1560	tip=1;
1561	oldOlist=list(0);
1562	oldC=list(0);
1563	oldt=list(0);
1564	oldD=list(0);
1565	oldH=list(0);
1566	allH=list(0); // the value of the parent + the new one
1567	infobo7=-1;
1568
1569	auxmb=tip,oldOlist,oldC,oldt,oldD,oldH,allH,infobo7;
1570	mobile[idchart+1]=auxmb;
1571
1572	idchart=idchart+1;}
1573
1574	else{ // CASE EordI==0 AND I!=1 --> HYPERBOLIC EQUATIONS
1575
1576	// COMPUTE THE IDEAL OF NON HYPERBOLIC ELEMENTS
1577
1578	extraL=Nonhyp(Coef,expI,k,n,flag,auxordI[2]); // gives (newcoef,newI,hypcoef,genhyp,flaglist)
1579
1580	// CHECK IF ALL THE VARIABLES ARE ALREADY INVERTIBLE
1581
1582	newflag=extraL[5];
1583	chart[idchart][6]=extraL[5]; // update the status of variables
1584
1585	cont3=0;
1586	for (i=1;i<=n;i++){if (newflag[i]==1){cont3=cont3+1;}}
1587
1588	if (cont3==n){ // ALL THE VARIABLES ARE INVERTIBLE
1589	auxfchart[1]=chart[idchart];
1590	finalchart=finalchart+auxfchart;
1591	empstring="@";
1592	invariant[idchart]=empstring;
1593	auxinv[idchart]=list(0);
1594	nsons[idchart]=0;
1595
1596	// information for the mobile
1597	tip=1;
1598	oldOlist=list(0);
1599	oldC=list(0);
1600	oldt=list(0);
1601	oldD=list(0);
1602	oldH=list(0);
1603	allH=list(0);
1604	infobo7=-1;
1605
1606	auxmb=tip,oldOlist,oldC,oldt,oldD,oldH,allH,infobo7;
1607	mobile[idchart+1]=auxmb;
1608
1609	idchart=idchart+1;}
1610	else{ // OTHERWISE, CONTINUE CHEKING IF newI=0 or not
1611
1612	Coef=extraL[1];
1613	expI=extraL[2];
1614
1615	hipercoefson=extraL[3]; // Information about hyperbolic generators
1616	hiperexpson=extraL[4];
1617
1618	k=size(expI);
1619
1620	if (k==0){auxfchart[1]=chart[idchart]; // WE HAVE FINISHED
1621	finalchart=finalchart+auxfchart;
1622	empstring="#"; // no more non-hyperbolic generators in this chart
1623	invariant[idchart]=empstring;
1624	auxinv[idchart]=list(0);
1625	nsons[idchart]=0;
1626
1627	// information for the mobile
1628	tip=1;
1629	oldOlist=list(0);
1630	oldC=list(0);
1631	oldt=list(0);
1632	oldD=list(0);
1633	oldH=list(0);
1634	allH=list(0);
1635	infobo7=-1;
1636
1637	auxmb=tip,oldOlist,oldC,oldt,oldD,oldH,allH,infobo7;
1638	mobile[idchart+1]=auxmb;
1639
1640	idchart=idchart+1;}
1641
1642	else{ // CONTINUE WITH THE IDEAL OF NON HYPERBOLIC EQS
1643
1644	chart[idchart][4]=expI; // new input ideal and coefficients
1645	chart[idchart][5]=Coef;
1646	chart[idchart][10]=hipercoef+hipercoefson;
1647	chart[idchart][11]=hiperexp+hiperexpson;
1648
1649	idchart=idchart;
1650	control2=1; // it is necessary to reset the value of c
1651	control3=1; // and the previous exceptional divisors
1652	}
1653
1654	// PROBABLY IT IS NEC MORE INFORMATION !!!
1655
1656	} // closing else otherwise
1657
1658	} // closing else case I!=1
1659
1660	} // closing else for EordI=0
1661
1662	if (EordI<0) {print("ERROR in chart"); print(idchart); break;}
1663
1664	//----------------------- guardar de momento--------
1665	// if (EordI==0) {auxfchart[1]=chart[idchart];
1666	// finalchart=finalchart+auxfchart;
1667	// L2=Gamma(expJ,c,n); // HAY QUE APLICARLO AL M NO AL J
1668	// invariant[idchart]=L2[2];
1669	// auxinv[idchart]=list(0);
1670	// nsons[idchart]=0;
1671	// idchart=idchart+1;}
1672	//------------------------------------------------
1673
1674
1675	} // END ELSE
1676	//---------------------------------------------------
1677
1678	} // END LOOP WHILE
1679
1680	tnum=step[nchart];
1681
1682	totalinv=resfunction(invariant,auxinv,nchart,n);
1683
1684	return(chart,finalchart,invariant,nchart,step,nsons,auxinv,mobile,totalinv);
1685	}
1686	example
1687	{"EXAMPLE:"; echo = 2;
1688	ring r = 0,(x(1..2)),dp;
1689	ideal J=x(1)^2-x(2)^3;
1690	list L=Eresol(J);
1691	L[1][1]; // information of the first chart, L[1] list of charts
1692	L[2]; // list of charts with information about sons
1693	L[3]; // invariant, "#" means solved chart
1694	L[4]; // number of charts, 7 in this example
1695	L[5]; // height corresponding to each chart
1696	L[6]; // number of sons
1697	L[7]; // auxiliary invariant
1698	L[8]; // H exceptional divisors and more information
1699	L[9]; // complete resolution function
1700
1701	"Second example, write L[i] to see the i-th component of the list";
1702	ring r = 0,(x(1..3)),dp;
1703	ideal J=x(1)^2*x(2),x(3)^3; // SOLVED!
1704	list L=Eresol(J);
1705	L[4]; // 16 charts
1706	L[9]; // complete resolution function
1707
1708	"Third example, write L[i] to see the i-th component of the list";
1709	ring r = 0,(x(1..2)),dp;
1710	ideal J=x(1)^3-x(1)*x(2)^3;
1711	list L=Eresol(J);
1712	L[4]; // 8 charts, rational exponents
1713	L[9]; // complete resolution function
1714	}
1715
1716	//////////////////////////////////////////////////////////////////////////////////////
1717
1718	proc resfunction(list invariant, list auxinv, int nchart,int n)
1719	"USAGE: resfunction(invariant,auxinv,nchart,n);
1720	invariant, auxinv lists, nchart, n integers
1721	COMPUTE: Patch the resolution function
1722	RETURN: The complete resolution function
1723	EXAMPLE: example resfunction; shows an example
1724	"
1725	{
1726	int i,j,l,k;
1727	list patchfun,aux;
1728
1729	for (i=1;i<=nchart;i++){patchfun[i]=invariant[i];}
1730
1731	for (i=1;i<=nchart;i++){if (auxinv[i][1]!=0 and size(auxinv[i])==3){l=size(invariant[i]);
1732	for (j=1;j<=l;j++){
1733	if (invariant[i][j]==0){aux=auxinv[i];
1734	patchfun[i][j]=aux;
1735	if (l<n){for (k=j+1;k<=n;k++){patchfun[i][k]="*";}}}}
1736
1737	}
1738	else{
1739	if (auxinv[i][1]==1 and size(auxinv[i])==1){l=size(invariant[i]);
1740	if (l<n){for (k=l+1;k<=n;k++){patchfun[i][k]="*";}}
1741	}
1742	}
1743	}
1744
1745	return(patchfun);
1746	}
1747	example
1748	{"EXAMPLE:"; echo = 2;
1749	ring r = 0,(x(1..2)),dp;
1750	ideal J=x(1)^2-x(2)^3;
1751	list L=Eresol(J);
1752	L[3]; // incomplete resolution function
1753	resfunction(L[3],L[7],7,2); // complete resolution function
1754	}
1755	//////////////////////////////////////////////////////////////////////////////////////
1756	// //
1757	// MAIN PROCEDURE //
1758	// //
1759	//////////////////////////////////////////////////////////////////////////////////////
1760
1761	proc BINresol(ideal J)
1762	"USAGE: BINresol(J); J ideal
1763	RETURN: E-resolution of singularities of a binomial ideal J in terms of the affine charts, see example
1764	EXAMPLE: example BINresol; shows an example
1765	"
1766	{
1767
1768	int p,n;
1769
1770	p=char(basering);
1771	n=nvars(basering); // YA SE CALCULA EN Eresol, MEJORAR?
1772
1773	if (p>0){list Lring=ringlist(basering);
1774	Lring[1]=0;
1775	def r=basering;
1776	def Rnew=ring(Lring);
1777	setring Rnew;
1778	ideal chy=maxideal(1);
1779	map fRnew=r,chy;
1780	ideal J=fRnew(J);
1781
1782	// E-RESOLUTION, Computations in char 0
1783
1784	list L=Eresol(J);
1785
1786	// STEP 2: WRITE THE LOCALLY MONOMIAL IDEAL AS A MONOMIAL IDEAL
1787
1788	// not implemented yet, CHAR p !!!!
1789
1790	// STEP 3: DO THE E-RESOLUTION AGAIN (char 0 again)
1791
1792
1793	// generating output
1794
1795	int q=lcmofall(L[4],L[8]); // lcm of the denominators
1796
1797	list B=genoutput(L[1],L[8],L[4],L[6],n,q); // generate output needed for visualization
1798
1799
1800	// setring r; // Back to the basering
1801	// ideal chy=maxideal(1);
1802	// map fr=Rnew,chy;
1803	// list L=fr(L);
1804	// list B=fr(B);
1805
1806	}
1807
1808	else{
1809
1810	// E-RESOLUTION
1811
1812	list L=Eresol(J);
1813
1814	// STEP 2: WRITE THE LOCALLY MONOMIAL IDEAL AS A MONOMIAL IDEAL
1815
1816	// not implemented yet
1817
1818	// STEP 3: DO THE E-RESOLUTION AGAIN
1819
1820
1821	// generating output
1822
1823	int q=lcmofall(L[4],L[8]);
1824
1825	list B=genoutput(L[1],L[8],L[4],L[6],n,q);
1826
1827	}
1828
1829	return(B);
1830	}
1831	example
1832	{"EXAMPLE:"; echo = 2;
1833	ring r = 0,(x(1..2)),dp;
1834	ideal J=x(1)^2-x(2)^3;
1835	list B=BINresol(J);
1836	B[1]; // list of final charts
1837	B[2]; // list of all charts
1838
1839	ring r2 = 2,(x(1..3)),dp;
1840	ideal J=x(1)^2-x(2)^2*x(3)^2;
1841	list B2=BINresol(J);
1842	B2[2]; // list of all charts
1843	}
1844	///////////////////////////////////////////////////////
1845
1846	proc Maxord(list L,int n)
1847	"USAGE: Maxord(L,n); L list, n integer
1848	COMPUTE: Find the maximal entry of a list, input is a list defining a monomial
1849	RETURN: maximum entry of a list and its position
1850	EXAMPLE: example Maxord; shows an example
1851	"
1852	{int i,can;
1853	number canmax;
1854	list aux;
1855
1856	canmax=1;
1857	can=1;
1858	for (i=1;i<=n;i++)
1859	{ if (L[i]>=canmax and i>=can)
1860	{canmax=L[i]; can=i;}}
1861
1862	return(canmax,can);
1863	}
1864	example
1865	{"EXAMPLE:"; echo = 2;
1866	ring r = 0,(x(1..3)),dp;
1867	ideal J=x(1)^2x(2)x(3)^5;
1868	list L=data(J,1,3);
1869	L[2]; // list of exponents
1870	Maxord(L[2][1][1],3);
1871	}
1872	///////////////////////////////////////////////////////
1873
1874	proc Gamma(list expM,number c,int n)
1875	"USAGE: Gamma(L,c,n); L list, c number, n integer
1876	COMPUTE: The Gamma function, resolution function corresponding to the monomial case
1877	RETURN: lists of maximum exponents in L, value of Gamma function, center of blow up
1878	EXAMPLE: example Gamma; shows an example
1879	"
1880	{int i,j,k,l,cont,can;
1881	intvec upla;
1882	number canmax;
1883	list expM2,gamma,L,aux,maxlist,center,aux2;
1884
1885	i=1;
1886	cont=0;
1887	expM2=expM;
1888
1889	while (cont==0 and i<=n)
1890	{
1891
1892	L=Maxord(expM2,n);
1893	aux=L[1];
1894	maxlist=maxlist + aux;
1895	can=L[2];
1896
1897	if (i==1) {upla=can; center=can;}
1898	else {upla=upla,can; aux2=can; center=center+aux2;}
1899
1900	canmax=sum(maxlist);
1901	if (canmax>=c)
1902	{gamma[1]=-i; gamma[2]=canmax/c; gamma[3]=upla; cont=1;}
1903	else {expM2[can]=0;}
1904	i=i+1;
1905	}
1906	return(maxlist,gamma,center);
1907	}
1908	example
1909	{"EXAMPLE:"; echo = 2;
1910	ring r = 0,(x(1..5)),dp;
1911	ideal J=x(1)^2x(2)x(3)^5x(4)^2x(5)^3;
1912	list L=data(J,1,5);
1913	list G=Gamma(L[2][1][1],9,5); // critical value c=9
1914	G[1]; // maximum exponents in the ideal
1915	G[2]; // maximal value of Gamma function
1916	G[3]; // center given by Gamma
1917	}
1918	///////////////////////////////////////////////////////
1919
1920	proc convertdata(list C,list L, int n, list flaglist)
1921	"USAGE: convertdata(C,L,n,flaglist);
1922	C, L, flaglist lists, n integer
1923	COMPUTE: Compute the ideal corresponding to the given lists C,L
1924	RETURN: an ideal whose coefficients are given by C, exponents given by L
1925	EXAMPLE: example convertdata; shows an example
1926	"
1927	{int i,j,k,a,b,lon;
1928	poly aux,aux1,aux2,aux3,f;
1929	ideal J;
1930
1931	aux=poly(0);
1932	aux1=poly(1);
1933	aux2=poly(0);
1934	aux3=poly(1);
1935
1936
1937	k=size(L);
1938	for (i=1;i<=k;i++){lon=size(C[i]);
1939	if (lon==1){ // variables in the monomial
1940	for (j=1;j<=n;j++){a=int(poly(L[i][1][j]));
1941	if (a!=0){
1942	if (flaglist[j]==0){aux=poly(x(j)^a);
1943	aux1=aux1*aux;}
1944	else {aux=poly(y(j)^a);
1945	aux1=aux1*aux;}
1946	}
1947	}
1948	if (C[i][1]!=0){aux1=C[i][1]*aux1;} // we add the coefficient
1949	else {aux1=0;}
1950
1951	J[i]=aux1;
1952	aux1=poly(1);
1953	}
1954
1955	else{ // variables in the binomial
1956
1957	for (j=1;j<=n;j++){a=int(poly(L[i][1][j])); b=int(poly(L[i][2][j]));
1958
1959	if (a!=0){
1960	if (flaglist[j]==0){aux=poly(x(j)^a);
1961	aux1=aux1*aux;}
1962	else {aux=poly(y(j)^a);
1963	aux1=aux1*aux;}
1964	}
1965
1966	if (b!=0){
1967	if (flaglist[j]==0){aux2=poly(x(j)^b);
1968	aux3=aux3*aux2;}
1969	else {aux2=poly(y(j)^b);
1970	aux3=aux3*aux2;}
1971	}
1972	}
1973	// we add the coefficients
1974	if (C[i][1]!=0){aux1=C[i][1]*aux1;}
1975	else {aux1=0;}
1976	if (C[i][2]!=0){aux3=C[i][2]*aux3;}
1977	else {aux3=0;}
1978
1979	f=aux1+aux3;
1980	J[i]=f;
1981	aux1=poly(1);
1982	aux3=poly(1);
1983
1984	}
1985	}
1986	return(J);
1987	}
1988	example
1989	{"EXAMPLE:"; echo = 2;
1990	ring r = 0,(x(1..4),y(5)),dp;
1991	list M=identifyvar();
1992	ideal J=x(1)^2y(5)^2-x(2)^2x(3)^2,6*x(4)^2;
1993	list L=data(J,2,5);
1994	L[1]; // Coefficients
1995	L[2]; // Exponents
1996	ideal J2=convertdata(L[1],L[2],5,M);
1997	J2;
1998	}
1999
2000	/////////////////////////////////////////////////////////////////////////////
2001
2002	proc lcmofall(int nchart,list mobile)
2003	"USAGE: lcmofall(nchart,mobile);
2004	nchart integer, mobile list of lists
2005	COMPUTE: Compute the lcm of the denominators of the E-orders of all the charts
2006	RETURN: an integer given the lcm
2007	NOTE: CALL BEFORE salida
2008	EXAMPLE: example lcmofall; shows an example
2009	"
2010
2011	{
2012	int i,m,tip,mcmall;
2013	intvec numall;
2014
2015	for (i=2;i<=nchart+1;i++){
2016	tip=mobile[i][1];
2017	if (tip!=1){numall=numall,tip;}
2018	}
2019	m=size(numall);
2020
2021	if (m==1){mcmall=1;}
2022	else{
2023	if (numall[1]==0){numall=numall[2..m];}
2024	mcmall=lcm(numall);}
2025
2026	return(mcmall);
2027	}
2028	example
2029	{"EXAMPLE:"; echo = 2;
2030	ring r = 0,(x(1..2)),dp;
2031	ideal J=x(1)^3-x(1)*x(2)^3;
2032	list L=Eresol(J);
2033	L[4]; // 8 charts, rational exponents
2034	L[8][2][2]; // E-orders at the first chart
2035	lcmofall(8,L[8]);
2036	}
2037	/////////////////////////////////////////////////////////////////////////////
2038
2039	proc salida(int idchart,list chart,list mobile,int numson,intvec previousa,int n,int q)
2040	"USAGE: salida(idchart,chart,mobile,numson,previousa,n,q);
2041	idchart, numson, n, q integers, chart, mobile, lists, previousa intvec
2042	COMPUTE: CONVERT THE OUTPUT OF A CHART IN A RING, WHERE DEFINE A BASIC OBJECT (BO)
2043	RETURN: the ring corresponding to the chart
2044	EXAMPLE: example salida; shows an example
2045	"
2046	{
2047	int l,i,m,aux,parent,m4,j;
2048	intvec Hhist,EOhist,aux7,aux9;
2049	list expJ,Coef,BO,blwhist,Eolist,hipercoef,hiperexp;
2050	list flag;
2051
2052	// chart gives: parent,Y,a,expJ,Coef,flag,Hhist,blwhist,path,hipercoef,hiperexp
2053	// mobile gives: tip,oldOlist,oldC,oldt,oldD,oldH,allH,infobo7; NOTE: Eolist=mobile[2];
2054
2055	// we need to define the suitable ring at this chart
2056
2057	list Lring=ringlist(basering);
2058	def RR2=basering;
2059
2060	flag=chart[6];
2061	string newl;
2062
2063	for (l=1;l<=n; l++){if (flag[l]==1){newl=string(l);
2064	Lring[2][l]="y("+newl+")";} }
2065
2066
2067	def RRnew=ring(Lring);
2068	setring RRnew;
2069	ideal chy=maxideal(1);
2070	map fRnew=RR2,chy;
2071
2072	list chart=fRnew(chart);
2073
2074	list mobile2=fRnew(mobile);
2075
2076
2077	flag=chart[6];
2078
2079	// we need to convert expJ and Coef to an ideal
2080
2081	expJ=chart[4];
2082	Coef=chart[5];
2083	Hhist=chart[7];
2084	blwhist=chart[8];
2085
2086	// now the ideal will be correctly defined in the ring Rnew
2087
2088	ideal J2=convertdata(Coef,expJ,n,flag); // Computations in RRnew
2089
2090	//------------------------------------------------------------------------------
2091	// START TO CREATE THE BO corresponding to this chart
2092
2093	BO=createBO(J2);
2094
2095	// MODIFY BO WITH THE INFORMATION OF THE CHART
2096
2097	// BO[1] an ideal, say W_i, defining the ambient space of the i-th chart of the blowing up
2098	// If there are hyperbolic equations, we put them here
2099
2100	hipercoef=chart[10];
2101	hiperexp=chart[11];
2102
2103	if (size(hipercoef)!=0){
2104	ideal ambJ=convertdata(hipercoef,hiperexp,n,flag);
2105	BO[1]=ambJ;
2106	}
2107
2108	// BO[2] an ideal defining the controlled transform
2109
2110	BO[2]=J2;
2111
2112	// BO[3] intvec, tupla containing the maximal E-order of BO[2]
2113
2114	if (numson==0){BO[3]=1;} // we write 1 if the chart is a final chart
2115	else{
2116	Eolist=mobile2[2]; // otherwise, convert the list of E-orders in an intvec
2117	m=size(Eolist);
2118	aux=int(Eolist[1]*q);
2119	EOhist=aux;
2120
2121	if (m>1){for (i=2;i<=m;i++){aux=int(Eolist[i]*q); EOhist=EOhist,aux;}}
2122
2123	BO[3]=EOhist;
2124	}
2125
2126	// BO[4] the list of exceptional divisors given by Hhist
2127
2128	BO[4]=constructH(Hhist,n,flag);
2129
2130	// BO[5] an ideal defining the map K[W] ----> K[Wi] given by blwhist
2131
2132	BO[5]=constructblwup(blwhist,n,chy,flag);
2133
2134	// BO[6] an intvec, BO[6][j]=1 indicates that <BO[4][j],BO[2]>=1, i.e. the
2135	// strict transform does not meet the j-th exceptional divisor
2136
2137	m4=size(BO[4]);
2138	ideal auxydeal;
2139	ideal Jint;
2140
2141	for (j=1;j<=m4;j++){
2142
2143	auxydeal=BO[4][j]+J2;
2144	Jint=std(auxydeal);
2145
2146	if (size(Jint)==1 and Jint[1]==1){BO[6][j]=1;}
2147	else{BO[6][j]=0;}
2148	}
2149
2150	// BO[7] intvec, the index of the first blown-up object in the resolution process
2151	// leading to this object for which the value of b was BO[3]
2152	// the subsequent ones are the indices for the Coeff-Objects
2153	// of BO[2] used when determining the center
2154	// index of last element of H^- in H
2155
2156
2157	if (numson!=0){BO[7]=mobile2[8];} // it is always -1 at the final charts
2158
2159	// BO[8] a matrix indicating that BO[4][i] meets BO[4][j] by BO[8][i,j]=1 for i < j
2160
2161	if (m4>0){
2162	matrix aux8[m4][m4];
2163
2164	BO[8]=aux8;
2165
2166	ideal auxydeal2;
2167	ideal Jint2;
2168
2169	for (i=1;i<=m4;i++){
2170	for (j=i+1;j<=m4;j++){
2171	auxydeal2=BO[4][i]+BO[4][j];
2172	Jint2=std(auxydeal2);
2173
2174	if (size(Jint2)==1 and Jint2[1]==1){BO[8][i,j]=0;}
2175	else{ for (l=1;l<j;l++){BO[8][l,j]=1;} }
2176	}
2177
2178	}
2179	}
2180	else{ matrix aux8[1][1];
2181	BO[8]=aux8;}
2182
2183
2184	// BO[9] INTERNAL DATA, second component of Villamayor resolution function,
2185	// only needed to use the visualization procedures
2186
2187	int m3=size(BO[3]);
2188
2189	if (m3==1){aux9=intvec(0);}
2190	else{ aux9[1]=0;
2191	for (i=2;i<=m3;i++){aux9=aux9,0;}
2192	}
2193
2194	BO[9]=aux9;
2195
2196	//------------------------------------------------------------------------------
2197
2198	// START TO CREATE THE extra information corresponding to this chart
2199
2200	/////////////// Short description of data in a chart ///////////////////
2201	// All chart data is stored in an object of type ring, the following
2202	// variables are always present in such a ring:
2203
2204	// BO: already created
2205
2206	// cent: ideal, describing the upcoming center determined by the algorithm
2207
2208	ideal cent=tradtoideal(previousa,J2,flag);
2209	export cent;
2210
2211	// path= module (autoconverted to matrix)
2212	// path[1][idchart]=parent[idchart] index of the parent-chart in resolution history of this chart
2213	// path[2][idchart]=index of this chart in relation with its brother-charts
2214
2215	module path=chart[9];
2216	export path;
2217
2218	// lastMap: ideal, describing the preceding blow up leading to this chart
2219
2220	ideal lastMap=constructlastblwup(blwhist,n,chy,flag);
2221	export lastMap;
2222
2223	//------------------------------------------------------------------------------
2224
2225	// EXTRA INFORMATION NEEDED
2226
2227	list invSat=ideal(0),aux9;
2228	export(invSat);
2229
2230	export BO;
2231
2232	return(RRnew);
2233	}
2234	example
2235	{"EXAMPLE:"; echo = 2;
2236	ring r = 0,(x(1..2)),dp;
2237	ideal J=x(1)^2-x(2)^3;
2238	list L=Eresol(J);
2239	list B=salida(5,L[1][5],L[8][6],2,L[1][3][3],2,1); // chart 5
2240	def RR=B[1];
2241	setring RR;
2242	BO;
2243
2244	ring r = 0,(x(1..2)),dp;
2245	ideal J=x(1)^2-x(2)^3;
2246	list L=Eresol(J);
2247	list B=salida(7,L[1][7],L[8][8],0,L[1][5][3],2,1); // chart 7
2248	def RR=B[1];
2249	setring RR;
2250	BO;
2251	showBO(BO);
2252
2253	ring r = 0,(x(1..2)),dp;
2254	ideal J=x(1)^3-x(1)*x(2)^3;
2255	list L=Eresol(J); // 8 charts, rational exponents
2256	list B=salida(1,L[1][1],L[8][2],2,0,2,2); // CHART 1
2257	def RR=B[1];
2258	setring RR;
2259	BO;
2260
2261	}
2262
2263	/////////////////////////////////////////////////////////////////////////////
2264	// CONVERT THE OUTPUT OF Eresol IN A LIST OF RINGS, WHERE A BASIC OBJECT (BO) IS DEFINED
2265	// IN ORDER TO INTEGRATE THIS LIBRARY INSIDE THE LIBRARY resolve.lib
2266
2267	proc genoutput(list chart,list mobile,int nchart,list nsons,int n,int q)
2268	"USAGE: genoutput(chart,mobile,nchart,nsons,n,q);
2269	chart, mobile, nsons lists, nchart, n, q integers
2270	RETURN: two lists, the first one gives the rings corresponding to the final charts,
2271	the second one is the list of all rings corresponding to the affine charts of the resolution process
2272	EXAMPLE: example genoutput; shows an example
2273	"
2274	{
2275	int idchart,parent;
2276	list auxlist,solvedrings,totalringlist,previousa;
2277
2278	// chart gives: parent,Y,a,expJ,Coef,flag,Hhist,blwhist,path,hipercoef,hiperexp
2279	// mobile gives: tip,oldOlist,oldC,oldt,oldD,oldH,allH,infobo7; NOTE: Eolist=mobile[2];
2280
2281	idchart=1;
2282
2283	// first loop, construct list previousa
2284
2285	while (idchart<=nchart)
2286	{
2287	if (idchart==1){previousa[1]=chart[2][3];}
2288	else
2289	{
2290	// if there are no sons, the next center is nothing
2291	if (nsons[idchart]==0){previousa[idchart]=0;}
2292	// always fill the parent
2293	parent=chart[idchart][1];
2294	previousa[parent]=chart[idchart][3];
2295	}
2296	idchart=idchart+1;
2297	}
2298
2299	// HERE BEGIN THE LOOP
2300
2301	idchart=1;
2302
2303	while (idchart<=nchart)
2304	{
2305	def auxexit=salida(idchart,chart[idchart],mobile[idchart+1],nsons[idchart],previousa[idchart],n,q);
2306	// we add the ring to the list of all rings
2307	auxlist[1]=auxexit;
2308	totalringlist=totalringlist+auxlist;
2309	// if the chart has no sons, add it to the list of final charts
2310	if (nsons[idchart]==0){solvedrings=solvedrings+auxlist;}
2311	auxlist=list();
2312	kill auxexit;
2313	idchart=idchart+1;
2314	} // EXIT WHILE
2315	return(solvedrings,totalringlist);
2316	}
2317	example
2318	{"EXAMPLE:"; echo = 2;
2319	ring r = 0,(x(1..2)),dp;
2320	ideal J=x(1)^3-x(1)*x(2)^3;
2321	list L=Eresol(J); // 8 charts, rational exponents
2322	list B=genoutput(L[1],L[8],L[4],L[6],2,2); // generates the output
2323	presentTree(B);
2324	list iden0=collectDiv(B);
2325	ResTree(B,iden0[1]); // generates the resolution tree
2326
2327	// Use presentTree(B); to see the final charts
2328	// To see the tree type in another shell
2329	// dot -Tjpg ResTree.dot -o ResTree.jpg
2330	// /usr/bin/X11/xv ResTree.jpg
2331	}
2332	/////////////////////////////////////////////////////////////////////
2333
2334	proc computemcm(list Eolist)
2335	"USAGE: computemcm(Eolist); Eolist list
2336	RETURN: an integer, the least common multiple of the denominators of the E-orders
2337	NOTE: Make the same as lcmofall but for one chart. NECESSARY BECAUSE THE E-ORDERS ARE OF TYPE NUMBER!!
2338	EXAMPLE: example computemcm; shows an example
2339	"
2340	{
2341	int m,i,aux,mcmchart;
2342	intvec num;
2343
2344	m=size(Eolist);
2345
2346	if (m==1)
2347	{
2348	mcmchart=int(denominator(Eolist[1]));
2349	return(mcmchart);
2350	}
2351
2352	if (m>1)
2353	{
2354	num=int(denominator(Eolist[1]));
2355	for (i=2;i<=m;i++)
2356	{
2357	aux=int(denominator(Eolist[i]));
2358	num=num,aux;
2359	}
2360	}
2361	mcmchart=lcm(num);
2362	return(mcmchart);
2363	}
2364	example
2365	{"EXAMPLE:"; echo = 2;
2366	ring r = 0,(x(1..2)),dp;
2367	ideal J=x(1)^3-x(1)*x(2)^3;
2368	list L=Eresol(J); // 8 charts, rational exponents
2369	L[8][2][2]; // maximal E-order at the first chart
2370	computemcm(L[8][2][2]);
2371	}
2372	/////////////////////////////////////////////////////////////////////
2373
2374	proc constructH(intvec Hhist,int n,list flag)
2375	"USAGE: constructH(Hhist,n,flag);
2376	Hhist intvec, n integer, flag list
2377	RETURN: the list of exceptional divisors accumulated at this chart
2378	EXAMPLE: example constructH; shows an example
2379	"
2380	{
2381	int i,j,m,l;
2382	list exceplist;
2383	ideal aux;
2384
2385	m=size(Hhist);
2386	if (Hhist[1]==0 and m>1){Hhist=Hhist[2..m]; m=m-1;
2387
2388	for (i=1;i<=m;i++){
2389	l=Hhist[i];
2390	if (flag[l]==0){aux=ideal(poly(x(l))); }
2391	else {aux=ideal(poly(y(l))); }
2392
2393	exceplist[i]=aux;
2394	}
2395	// eliminate repeated variables
2396	for (i=1;i<=m;i++){for (j=1;j<=m;j++){
2397	if (Hhist[i]==Hhist[j] and i!=j){
2398	if (i<j){exceplist[i]=ideal(1);}
2399	if (i>j){exceplist[j]=ideal(1);}
2400	}
2401	}
2402	}
2403
2404	}
2405	else {exceplist=list();}
2406
2407	// else {exceplist=list(ideal(0));} // IF IT FAILS USE THIS
2408
2409	return(exceplist);
2410	}
2411	example
2412	{"EXAMPLE:"; echo = 2;
2413	ring r = 0,(x(1..3)),dp;
2414	list flag=identifyvar();
2415	ideal J=x(1)^4*x(2)^2, x(1)^2+x(3)^3;
2416	list L=Eresol(J); // 7 charts
2417	// history of the exceptional divisors at the 7-th chart
2418	L[1][7][7]; // blow ups at x(3)-th, x(1)-th and x(1)-th charts
2419	constructH(L[1][7][7],3,flag);
2420
2421	}
2422	/////////////////////////////////////////////////////////////////////
2423
2424	proc constructblwup(list blwhist,int n,ideal chy,list flag)
2425	"USAGE: constructblwup(blwhist,n,chy,flag);
2426	blwhist, flag lists, n integer, chy ideal
2427	RETURN: the ideal defining the map K[W] --> K[Wi],
2428	which gives the composition map of all the blowing up leading to this chart
2429	NOTE: NECESSARY START WITH COLUMNS
2430	EXAMPLE: example constructblwup; shows an example
2431	"
2432	{
2433	int i,j,m,m2;
2434	poly aux2;
2435
2436	m=size(blwhist[1]);
2437
2438	for (j=1;j<=m;j++){
2439	for (i=1;i<=n;i++){ m2=blwhist[i][j];
2440
2441	// If m2!=0 this variable changes. First decide if the variable to multiply is invertible or not
2442
2443	if (m2!=0){
2444	if (flag[m2]==0){aux2=poly(x(m2));}
2445	else {aux2=poly(y(m2));}
2446
2447	// And then substitute this variable for the corresponding product in the whole ideal
2448
2449	if (flag[i]==0){chy=subst(chy,x(i),x(i)*aux2);}
2450	else {chy=subst(chy,y(i),y(i)*aux2);}
2451
2452	}
2453	}
2454
2455	}
2456
2457	return(chy);
2458	}
2459	example
2460	{"EXAMPLE:"; echo = 2;
2461	ring r = 0,(x(1..3)),dp;
2462	list flag=identifyvar();
2463	ideal chy=maxideal(1);
2464	ideal J=x(1)^4*x(2)^2, x(1)^2+x(3)^3;
2465	list L=Eresol(J); // 7 charts
2466	// history of the blow ups at the 7-th chart, center {x(1)=x(3)=0} every time
2467	L[1][7][8]; // blow ups at x(3)-th, x(1)-th and x(1)-th charts
2468	constructblwup(L[1][7][8],3,chy,flag);
2469	}
2470	/////////////////////////////////////////////////////////////////////
2471
2472	proc constructlastblwup(list blwhist,int n,ideal chy,list flag)
2473	"USAGE: constructlastblwup(blwhist,n,chy,flag);
2474	blwhist, flag lists, n integer, chy ideal
2475	RETURN: the ideal defining the last blow up
2476	NOTE: NECESSARY START WITH COLUMNS
2477	EXAMPLE: example constructlastblwup; shows an example
2478	"
2479	{
2480	int i,j,m,m2;
2481	poly aux2;
2482
2483	m=size(blwhist[1]);
2484
2485	if (m>0){
2486	for (i=1;i<=n;i++){ m2=blwhist[i][m];
2487
2488	// If m2!=0 this variable changes. First decide if the variable to multiply is invertible or not
2489
2490	if (m2!=0){
2491	if (flag[m2]==0){aux2=poly(x(m2));}
2492	else {aux2=poly(y(m2));}
2493
2494	// And then substitute this variable for the corresponding product in the whole ideal
2495
2496	if (flag[i]==0){chy=subst(chy,x(i),x(i)*aux2);}
2497	else {chy=subst(chy,y(i),y(i)*aux2);}
2498
2499	}
2500	}
2501	}
2502
2503	return(chy);
2504	}
2505	example
2506	{"EXAMPLE:"; echo = 2;
2507	ring r = 0,(x(1..3)),dp;
2508	list flag=identifyvar();
2509	ideal chy=maxideal(1);
2510	ideal J=x(1)^4*x(2)^2, x(1)^2+x(3)^3;
2511	list L=Eresol(J); // 7 charts
2512	// history of the blow ups at the 7-th chart, center {x(1)=x(3)=0} every time
2513	L[1][7][8]; // blow ups at x(3)-th, x(1)-th and x(1)-th charts
2514	constructlastblwup(L[1][7][8],3,chy,flag);
2515	}
2516	/////////////////////////////////////////////////////////////////////
2517
2518	proc tradtoideal(intvec a,ideal J2,list flag)
2519	"USAGE: tradtoideal(a,J2,flag);
2520	a intvec, J2 ideal, flag list
2521	COMPUTE: traslate to an ideal the intvec defining the center
2522	RETURN: the ideal of the center, given by the intvec a, or J2 if a=0
2523	EXAMPLE: example tradtoideal; shows an example
2524	"
2525	{
2526	int i,m;
2527	ideal acenter,aux2;
2528
2529	if (a==0){acenter=J2;}
2530	else{
2531	m=size(a);
2532	for (i=1;i<=m;i++){
2533	if (flag[a[i]]==0){aux2=poly(x(a[i]));}
2534	else {aux2=poly(y(a[i]));}
2535
2536	acenter=acenter+aux2;
2537	}
2538	}
2539	return(acenter);
2540	}
2541	example
2542	{"EXAMPLE:"; echo = 2;
2543	ring r = 0,(x(1..3)),dp;
2544	list flag=identifyvar();
2545	ideal J=x(1)^4*x(2)^2, x(1)^2+x(3)^3;
2546	intvec a=1,3; // first center of blowing up
2547	tradtoideal(a,J,flag);
2548	}
2549	//////////////////////////////////////////////////////////////////////////////////////
2550	// OPERATIONS WITH LISTS
2551	//////////////////////////////////////////////////////////////////////////////////////
2552
2553	proc iniD(int n)
2554	"USAGE: iniD(n); n integer
2555	RETURN: list of lists of zeros of size n
2556	EXAMPLE: example iniD; shows an example
2557	"
2558	{int i,j;
2559	list D,auxD;
2560	for (j=1;j<=n; j++) {auxD[j]=0;}
2561	for (i=1;i<=n; i++) {D[i]=auxD;}
2562	return(D);
2563	}
2564	example
2565	{"EXAMPLE:"; echo = 2;
2566	iniD(3);
2567	}
2568	/////////////////////////////////////////////////////////
2569
2570	proc sumlist(list L1,list L2)
2571	"USAGE: sumlist(L1,L2); L1,L2 lists, (size(L1)==size(L2))
2572	RETURN: a list, sum of L1 and L2
2573	EXAMPLE: example sumlist; shows an example
2574	"
2575	{
2576	int i,k;
2577	list sumL;
2578	k=size(L1);
2579	if (size(L2)!=k) {return("ERROR en sumlist, lists must have the same size");}
2580	for (i=1;i<=k;i++) {sumL[i]=L1[i]+L2[i];}
2581	return(sumL);
2582	}
2583	example
2584	{"EXAMPLE:"; echo = 2;
2585	list L1=1,2,3;
2586	list L2=5,9,7;
2587	sumlist(L1,L2);
2588	}
2589	///////////////////////////////////////////////////////
2590
2591	proc reslist(list L1,list L2)
2592	"USAGE: reslist(L1,L2); L1,L2 lists, (size(L1)==size(L2))
2593	RETURN: a list, subtraction of L1 and L2
2594	EXAMPLE: example reslist; shows an example
2595	"
2596	{
2597	int i,k;
2598	list resL;
2599	k=size(L1);
2600	if (size(L2)!=k) {return("ERROR en reslist, lists must have the same size");}
2601	for (i=1;i<=k;i++) {resL[i]=L1[i]-L2[i];}
2602	return(resL);
2603	}
2604	example
2605	{"EXAMPLE:"; echo = 2;
2606	list L1=1,2,3;
2607	list L2=5,9,7;
2608	reslist(L1,L2);
2609	}
2610	//////////////////////////////////////////////////////
2611
2612	proc multiplylist(list L,number a)
2613	"USAGE: multiplylist(L,a); L list, a number
2614	RETURN: list of elements of type number, multiplication of L times a
2615	EXAMPLE: example multiplylist; shows an example
2616	"
2617	{int i,k;
2618	list newL,bb;
2619	number b;
2620	k=size(L);
2621	for (i=1;i<=k;i++) {b=L[i]*a; bb=b; newL=newL+bb;}
2622	return(newL);
2623	}
2624	example
2625	{"EXAMPLE:"; echo = 2;
2626	ring r = 0,(x(1..3)),dp;
2627	list L=1,2,3;
2628	multiplylist(L,1/5);
2629	}
2630	///////////////////////////////////////////////////////
2631
2632	proc dividelist(list L1,list L2)
2633	"USAGE: dividelist(L1,L2); L1,L2 lists
2634	RETURN: list of elements of type number, division of L1 by L2
2635	EXAMPLE: example dividelist; shows an example
2636	"
2637	{int i,k,k1,k2;
2638	list LL,bb;
2639	number a1,a2,b;
2640	k1=size(L1);
2641	k2=size(L2);
2642	if (k2!=k1) {print("ERROR en dividelist, lists must have the same size");}
2643	if (k1<=k2) {k=k1;}
2644	else {k=k2;}
2645	for (i=1;i<=k;i++)
2646	{a1=L1[i]; a2=L2[i]; b=a1/a2; bb=b; LL=LL+bb;}
2647	return(LL);
2648	}
2649	example
2650	{"EXAMPLE:"; echo = 2;
2651	ring r = 0,(x(1..3)),dp;
2652	list L1=1,2,3;
2653	list L2=5,9,7;
2654	dividelist(L1,L2);
2655	}
2656	///////////////////////////////////////////////////////
2657
2658	proc createlist(list L1,list L2)
2659	"USAGE: createlist(L1,L2); L1,L2 lists, (size(L1)==size(L2))
2660	RETURN: list of lists of two elements, the first one of L1 and the second of L2
2661	EXAMPLE: example createlist; shows an example
2662	"
2663	{
2664	int i,k;
2665	list L,aux;
2666	k=size(L1);
2667	if (size(L2)!=k)
2668	{ERROR ("createlist: lists must have the same size");}
2669	L=list0(k);
2670	for (i=1;i<=k;i++)
2671	{
2672	if (L1[i]!=0)
2673	{
2674	aux=L1[i],L2[i]; L[i]=aux;
2675	}
2676	else {L=delete(L,i);}
2677	}
2678	return(L);
2679	}
2680	example
2681	{"EXAMPLE:"; echo = 2;
2682	list L1=1,2,3;
2683	list L2=5,9,7;
2684	createlist(L1,L2);
2685	}
2686	///////////////////////////////////////////////////////
2687	proc list0(int n)
2688	"USAGE: list0(n); n integer
2689	RETURN: list of n zeros
2690	EXAMPLE: example list0; shows an example
2691	"
2692	{
2693	int i;
2694	list L0;
2695	for (i=1;i<=n;i++) {L0[i]=0;}
2696	return(L0);
2697	}
2698	example
2699	{"EXAMPLE:"; echo = 2;
2700	list0(4);
2701	}
2702	////////////////////////////////////////////////////////////
2703
2704	proc Emaxcont(list Coef,list Exp,int k,int n,list flag)
2705	"USAGE: Emaxcont(Coef,Exp,k,n,flag);
2706	Coef,Exp,flag lists, k,n, integers
2707	Exp is a list of lists of exponents, k=size(Exp)
2708	COMPUTE: Identify ALL the variables of E-maximal contact
2709	RETURN: a list with the indexes of the variables of E-maximal contact
2710	EXAMPLE: example Emaxcont; shows an example
2711	"
2712	{
2713	int i,j,lon;
2714	number maxEo;
2715	list L,sums,bx,maxvar;
2716
2717	L=maxEord(Coef,Exp,k,n,flag);
2718
2719	maxEo=L[1];
2720	sums=L[2];
2721
2722	if (maxEo>0)
2723	{
2724	for (i=1;i<=k; i++)
2725	{
2726	lon=size(sums[i]);
2727	if (lon==2)
2728	{
2729	if (sums[i][1]==maxEo) // variables of the first term
2730	{
2731	for (j=1;j<=n; j++)
2732	{
2733	if(Exp[i][1][j]!=0 and flag[j]==0)
2734	{
2735	bx=j; maxvar=maxvar + bx;
2736	}
2737	}
2738	}
2739	if (sums[i][2]==maxEo) // variables of the second term
2740	{
2741	for (j=1;j<=n; j++)
2742	{
2743	if(Exp[i][2][j]!=0 and flag[j]==0){bx=j; maxvar=maxvar + bx;}
2744	}
2745	}
2746	}
2747	else
2748	{
2749	if (sums[i][1]==maxEo)
2750	{
2751	for (j=1;j<=n; j++)
2752	{
2753	if(Exp[i][1][j]!=0 and flag[j]==0)
2754	{
2755	bx=j; maxvar=maxvar + bx;
2756	}
2757	}
2758	}
2759	}
2760	}
2761	}
2762	else {maxvar=list();}
2763
2764	// eliminating repeated terms
2765	maxvar=elimrep(maxvar);
2766	// It is necessary to check if flag[j]==0 in order to avoid the selection of y variables
2767	return(maxEo,maxvar);
2768	}
2769	example
2770	{"EXAMPLE:"; echo = 2;
2771	ring r = 0,(x(1),y(2),x(3),y(4),x(5..7),y(8)),dp;
2772	list flag=identifyvar();
2773	ideal J=x(1)^3x(3)-y(2)y(4)^2,x(5)y(2)-x(7)y(4)^2,x(6)^2(1-y(4)y(8)^5),x(7)^4*y(8)^2;
2774	list L=data(J,4,8);
2775	list hyp=Emaxcont(L[1],L[2],4,8,flag);
2776	hyp[1]; // max E-order=0
2777	hyp[2]; // There are no hypersurfaces of E-maximal contact
2778
2779	ring r = 0,(x(1),y(2),x(3),y(4),x(5..7),y(8)),dp;
2780	list flag=identifyvar();
2781	ideal J=x(1)^3x(3)-y(2)y(4)^2x(3),x(5)y(2)-x(7)y(4)^2,x(6)^2(1-y(4)y(8)^5),x(7)^4y(8)^2;
2782	list L=data(J,4,8);
2783	list hyp=Emaxcont(L[1],L[2],4,8,flag);
2784	hyp[1]; // the E-order is 1
2785	hyp[2]; // {x(3)=0},{x(5)=0},{x(7)=0} are hypersurfaces of E-maximal contact
2786	}
2787	///////////////////////////////////////////////////////

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

Download in other formats:

Original Format