Context Navigation

resbin.lib @ 7c0817

spielwiese

Last change on this file since 7c0817 was 7c0817, checked in by Hans Schoenemann <hannes@…>, 14 years ago
new library resbin.lib git-svn-id: file:///usr/local/Singular/svn/trunk@12819 2c84dea3-7e68-4137-9b89-c4e89433aadc
Property mode set to `100644`
File size: 90.7 KB

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

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

Download in other formats:

Original Format