Context Navigation

binresol.lib @ c4c127

spielwiese

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

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

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

Download in other formats:

Original Format