Context Navigation

resbinomial.lib @ 19cb509

spielwiese

Last change on this file since 19cb509 was 19cb509, checked in by Hans Schoenemann <hannes@…>, 12 years ago
id git-svn-id: file:///usr/local/Singular/svn/trunk@14063 2c84dea3-7e68-4137-9b89-c4e89433aadc
Property mode set to `100644`
File size: 89.7 KB

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

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

Download in other formats:

Original Format