Context Navigation

binresol.lib @ 555878

spielwiese

Last change on this file since 555878 was 555878, checked in by Hans Schönemann <hannes@…>, 14 years ago
fixed docu, some format stuff git-svn-id: file:///usr/local/Singular/svn/trunk@12518 2c84dea3-7e68-4137-9b89-c4e89433aadc
Property mode set to `100644`
File size: 86.0 KB

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

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

Download in other formats:

Original Format