Context Navigation

resolve.lib @ 0f0b30

spielwiese

Last change on this file since 0f0b30 was 0f0b30, checked in by Hans Schoenemann <hannes@…>, 8 years ago
fix: some place of ideal==int
Property mode set to `100644`
File size: 134.7 KB

Line
1	////////////////////////////////////////////////////////////////////////////
2	version="version resolve.lib 4.0.0.0 Jun_2013 "; // $Id$
3	category="Algebraic Geometry";
4	info="
5	LIBRARY: resolve.lib Resolution of singularities (Desingularization)
6	Algorithm of Villamayor
7	AUTHORS: A. Fruehbis-Krueger, anne@mathematik.uni-kl.de,
8	@* G. Pfister, pfister@mathematik.uni-kl.de
9
10	REFERENCES:
11	[1] J.Kollar: Lectures on Resolution of Singularities, Princeton University Press (2007)@*
12	(contains large overview over various known methods for curves and surfaces as well as@*
13	a detailed description of the approach in the general case)@*
14	[2] A.Bravo, S.Encinas, O.Villamayor: A Simplified Proof of Desingularisation and@*
15	Applications, Rev. Math. Iberoamericana 21 (2005), 349-458@*
16	(description of the algorithmic proof of desingularization in characteristic zero
17	which underlies this implementation)@*
18	[3] A.Fruehbis-Krueger: Computational Aspects of Singularities, in J.-P. Brasselet,
19	J.Damon et al.: Singularities in Geometry and Topology, World Scientific
20	Publishing, 253--327 (2007)@*
21	(chapter 4 contains a detailed discussion on algorithmic desingularization and
22	efficiency aspects thereof)
23
24	PROCEDURES:
25	blowUp(J,C[,W,E]) computes the blowing up of the variety V(J) (considered
26	as embedded in V(W)) in the (smooth) center V(C),
27	blowUp2(J,C) computes the blowing up of the variety V(J) in the
28	(possibly singular) center V(C)
29	Center(J[,W,E]) computes 'Villamayor'-center for blow up
30	resolve(J) computes the desingularization of the variety V(J)
31
32	showBO(BO) prints the content of a BO in more human readable form
33	presentTree(L) prints the final charts in more human readable form
34	showDataTypes() prints help text for output data types
35
36	blowUpBO(BO,C) computes the blowing up of the variety V(BO[1]) in the
37	center V(C). BO is a list (basic object), C is an ideal
38	createBO(J,W,E) creates basic object from input data
39	CenterBO(BO) computes the center for the next blow-up of the
40	given basic object
41	Delta(BO) apply the Delta-operator of [Bravo,Encinas,Villamayor]
42	DeltaList(BO) list of results of Delta^0 to Delta^bmax
43	";
44	LIB "elim.lib";
45	LIB "primdec.lib";
46	LIB "presolve.lib";
47	LIB "linalg.lib";
48	LIB "sing.lib";
49	///////////////////////////////////////////////////////////////////////////////
50	// Tasks:
51	// 1) optimization of the local case
52	// 2) optimization in Coeff
53	// 3) change invariant to represent coeff=1 case
54	///////////////////////////////////////////////////////////////////////////////
55
56	proc showDataTypes()
57	"USAGE: showDataTypes();
58	RETURN: nothing, only pretty printing of extended version of help text
59	EXAMPLE: none
60	"
61	{
62	" ";
63	"//////////////// Short description of data type BO ///////////////////";
64	"BO[1] an ideal, say Wi, defining the ambient space of the i-th chart";
65	" of the blowing up";
66	"BO[2] an ideal defining the strict transform";
67	"BO[3] intvec, the first integer b such that in the original object";
68	" (Delta^b(BO[2]))==1";
69	" the subsequent integers have the same property for Coeff-Objects";
70	" of BO[2] used when determining the center";
71	"BO[4] the list of exceptional divisors";
72	"BO[5] an ideal defining the map K[W] ----> K[Wi]";
73	"BO[6] an intvec BO[6][j]=1 indicates that <BO[4][j],BO[2]>=1, i.e. the";
74	" strict transform does not meet the j-th exceptional divisor";
75	"BO[7] intvec,";
76	" the index of the first blown-up object in the resolution process";
77	" leading to this object for which the value of b was BO[3]";
78	" the subsequent ones are the indices for the Coeff-Objects";
79	" of BO[2] used when determining the center";
80	"BO[i], i>7: internal data";
81	" ";
82	pause();
83	" ";
84	"///////////// Short description of data in a chart ///////////////////";
85	"All chart data is stored in an object of type ring, the following ";
86	"variables are always present in such a ring:";
87	"BO: list of type basic object containing the variety, ambient space,";
88	" exceptional divisors and further data (see previous page)";
89	"cent: ideal, describing the upcoming center determined by the algorithm";
90	"path: module (autoconverted to matrix)";
91	" path[1][i]= (i-1)st chart in resolution history of this chart";
92	" path[2][i]= index of chart of the blow up leading to i-th chart";
93	" in resolution history";
94	"lastMap: ideal, describing the preceding blow up leading to this chart";
95	" ";
96	pause();
97	" ";
98	"///////////// Short description of type resolution data //////////////";
99	"list L containing two lists:";
100	"L[1]: list of rings, containing the final charts";
101	"L[2]: list of rings, containing all charts created in the resolution";
102	" process";
103	"The most convenient way to view these data is the procedure presentTree()";
104	"from this library. Alternatively, it can be digested using tools from";
105	"the libraries reszeta.lib and resgraph.lib";
106	"//////////////////////////////////////////////////////////////////////";
107	}
108	///////////////////////////////////////////////////////////////////////////////
109	proc createBO(ideal J,list #)
110	"USAGE: createBO(J[,W][,E]);
111	@* J,W = ideals
112	@* E = list
113	ASSUME: J = ideal containing W ( W = 0 if not specified)
114	@* E = list of smooth hypersurfaces (e.g. exceptional divisors)
115	RETURN: list BO representing a basic object :
116	BO[1] ideal W, if W has been specified; ideal(0) otherwise
117	BO[2] ideal J
118	BO[3] intvec
119	BO[4] the list E of exceptional divisors if specified;
120	empty list otherwise
121	BO[5] an ideal defining the identity map
122	BO[6] an intvec
123	BO[7] intvec
124	BO[8] a matrix
125	entries 3,5,6,7,8 are initialized appropriately for use of CenterBO
126	and blowUpBO
127	EXAMPLE: example createBO; shows an example
128	"
129	{
130	ideal W;
131	list E;
132	ideal abb=maxideal(1);
133	intvec v;
134	intvec bvec;
135	intvec w=-1;
136	matrix intE;
137	if(size(#)>0)
138	{
139	if(typeof(#[1])=="ideal")
140	{
141	W=#[1];
142	}
143	if(typeof(#[1])=="list")
144	{
145	E=#[1];
146	}
147	if(size(#)>1)
148	{
149	if((typeof(#[2])=="list") && (size(E)==0))
150	{
151	E=#[2];
152	}
153	if((typeof(#[2])=="ideal") && (size(W)==0))
154	{
155	W=#[2];
156	}
157	}
158	}
159	list BO=W,J,bvec,E,abb,v,w,intE;
160	return(BO);
161	}
162	example
163	{"EXAMPLE:";
164	echo = 2;
165	ring R=0,(x,y,z),dp;
166	ideal J=x2-y3;
167	createBO(J,ideal(z));
168	}
169	///////////////////////////////////////////////////////////////////////////////
170
171	proc blowUp(ideal J,ideal C,list #)
172	"USAGE: blowUp(J,C[,W][,E]);
173	W,J,C = ideals,
174	E = list
175	ASSUME: J = ideal containing W ( W = 0 if not specified)
176	@* C = ideal containing J
177	@* E = list of smooth hypersurfaces (e.g. exceptional divisors)
178	NOTE: W the ideal of the ambient space, C the ideal of the center of
179	the blowup and J the ideal of the variety
180	Important difference to blowUp2:
181	- the ambient space V(W) is blown up and V(J) transformed in it
182	@* - V(C) is assumed to be non-singular
183	COMPUTE: the blowing up of W in C, the exceptional locus, the strict
184	transform of J and the blowup map
185	RETURN: list, say l, of size at most size(C),
186
187	l[i] is the affine ring corresponding to the i-th chart
188	each l[i] contains the ideals
189	- aS, ideal of the blownup ambient space
190	- sT, ideal of the strict transform
191	- eD, ideal of the exceptional divisor
192	- bM, ideal corresponding to the blowup map
193
194	l[i] also contains a list BO, which can best be viewed with showBO(BO)
195	detailed information on the data type BO can be viewed via the
196	command showDataTypes();
197	EXAMPLE: example blowUp; shows an example
198	"
199	{
200	def S=basering;
201	ideal W;
202	list E;
203	ideal abb=maxideal(1);
204	intvec v;
205	intvec bvec;
206	intvec w=-1;
207	matrix intE;
208	if(size(#)>0)
209	{
210	if(typeof(#[1])=="ideal")
211	{
212	W=#[1];
213	}
214	if(typeof(#[1])=="list")
215	{
216	E=#[1];
217	}
218	if(size(#)>1)
219	{
220	if((typeof(#[2])=="list") && (size(E)==0))
221	{
222	E=#[2];
223	}
224	if((typeof(#[2])=="ideal") && (size(W)==0))
225	{
226	W=#[2];
227	}
228	}
229	}
230	list BO=W,J,bvec,E,abb,v,w,intE;
231	int locaT;
232	export locaT;
233	list blow=blowUpBO(BO,C,0);
234	kill locaT;
235	int i;
236	for(i=1;i<=size(blow);i++)
237	{
238	def Q=blow[i];
239	setring Q;
240	ideal aS=BO[1];
241	ideal sT=BO[2];
242	ideal eD=BO[4][size(BO[4])];
243	ideal bM=BO[5];
244	kill lastMap;
245	kill thisChart;
246	export(aS);
247	export(sT);
248	export(eD);
249	export(bM);
250	blow[i]=Q;
251	setring S;
252	kill Q;
253	}
254	return(blow);
255	}
256	example
257	{"EXAMPLE:";
258	echo = 2;
259	ring R=0,(x,y),dp;
260	ideal J=x2-y3;
261	ideal C=x,y;
262	list blow=blowUp(J,C);
263	def Q=blow[1];
264	setring Q;
265	aS;
266	sT;
267	eD;
268	bM;
269	}
270	///////////////////////////////////////////////////////////////////////////////
271
272	proc blowUp2(ideal J,ideal C)
273	"USAGE: blowUp2(J,C);
274	J,C = ideals,
275	ASSUME: C = ideal containing J
276	NOTE: C the ideal of the center of the blowup and J the ideal
277	of the variety
278	Important differences to blowUp:
279	- V(J) itself is blown up, not the ambient space
280	- C is not assumed to be non-singular
281	COMPUTE: the blowing up of J in C, the exceptional locus and the blow-up
282	map
283	RETURN: list, say l, of size at most size(C),
284
285	l[i] is the affine ring corresponding to the i-th chart
286	each l[i] contains the ideals
287	- Jnew, ideal of the blownup J
288	- eD, ideal of the new exceptional divisor
289	- bM, ideal corresponding to the blowup map
290
291	EXAMPLE: example blowUp2; shows an example
292	"
293	{
294	//----------------------------------------------------------------------------
295	// Initialization and basic sanity checks of the input data
296	//----------------------------------------------------------------------------
297	int i,j;
298	def S=basering;
299	def laM=maxideal(1);
300	//--- number of generators of C should be as small as possible
301	def mstdC=mstd(C);
302	if(ncols(mstdC[1])<ncols(mstdC[2]))
303	{
304	C=mstdC[1];
305	}
306	else
307	{
308	C=mstdC[2];
309	}
310	C=simplify(interred(C),2);
311	//--- the empty set is not a good center ;-)
312	if(deg(lead(C[1]))==0)
313	{
314	ERROR("Your chosen center was the empty set. Exiting.");
315	}
316	//--- V(C) should be a subset of V(J)
317	if(size(reduce(J,std(C)))>0)
318	{
319	ERROR("V(J) does not contain V(C). Exiting.");
320	}
321	//---------------------------------------------------------------------------
322	// To compute the blowing up, we need to consider
323	// (K[x_1,\dots,x_n]/J)[tC[1],...,tC[m]]
324	// which we want to represent as a quotient of K[x_1,...,x_n,y_1,..,y_m]
325	// by an ideal obtained by elimination of the extra variable t.
326	// (in the comments n=nvars(S) and m=ncols(C))
327	//---------------------------------------------------------------------------
328	//--- set up rings for the elimination
329	string r1="ring R=("+charstr(basering)+"),(x(1.."+string(nvars(S));
330	r1=r1+"),y(1.."+string(ncols(C))+")),dp;";
331	string r1t="ring Rt=("+charstr(basering)+"),(x(1.."+string(nvars(S));
332	r1t=r1t+"),y(1.."+string(ncols(C))+"),t),(dp,dp(1));";
333	execute(r1); // ring for describing the transforms
334	execute(r1t); // like r1, but with additional variable t
335	def J=fetch(S,J);
336	def C=fetch(S,C);
337	//--- we need to eliminate t from J,y(1)-tC[1],...y(m)-tC[m]
338	ideal elId=J;
339	for(i=1;i<=ncols(C);i++)
340	{
341	elId=elId,y(i)-t*C[i];
342	}
343	elId=eliminate(elId,t); // ideal describing the transform of J
344	setring R; // get rid of t
345	def elId=fetch(Rt,elId);
346	def E=fetch(S,C); // determine exceptional divisor
347	E=E+elId;
348	def laM0=fetch(S,laM); // the blowup map
349	//-----------------------------------------------------------------------------
350	// The result is now represented in an A^n \times P^{m-1}, hence
351	// involving n+m variables. For further computations we would like
352	// to pass to charts to keep the total number of variables as low as
353	// possible
354	//----------------------------------------------------------------------------
355	list resList;
356	ideal Jsub;
357	ideal Esub;
358	ideal laM;
359	ideal testId;
360	map phi;
361	list templist;
362	for(i=1;i<=nvars(R)-nvars(S);i++)
363	{
364	//--- first pass elId and E on to the i-th chart
365	Jsub=std(subst(elId,y(i),1));
366	if(deg(Jsub[1])==0)
367	{
368	//--- transform does not meet this chart ==> ignore it
369	i++;
370	continue;
371	}
372	Esub=std(subst(E,y(i),1));
373	//--- now get rid of unnecessary variables
374	//--- first by appropriate coordinate changes
375	templist=elimpart(Jsub);
376	if(size(templist[2])>0)
377	{
378	phi=R,templist[5];
379	Jsub=phi(Jsub);
380	Jsub=simplify(interred(Jsub),2);
381	Esub=phi(Esub);
382	Esub=simplify(interred(Esub),2);
383	laM=phi(laM0);
384	}
385	else
386	{
387	laM=laM0;
388	}
389	//--- then by dropping non-occuring variables
390	testId=Jsub,Esub,laM;
391	templist=findvars(testId);
392	if(size(templist[1])<nvars(R))
393	{
394	for(j=1;j<=size(templist[4]);j++)
395	{
396	laM=subst(laM,var(templist[4][j]),0);
397	}
398	execute("ring Rnew=("+charstr(basering)+"),("+string(templist[1])+"),dp;");
399	ideal Jnew=imap(R,Jsub);
400	ideal eD=imap(R,Esub);
401	ideal bM=imap(R,laM);
402	}
403	else
404	{
405	def Rnew=basering;
406	ideal Jnew=Jsub;
407	ideal eD=Esub;
408	ideal bM=laM;
409	}
410	//--- export the relevant data of this ring and add the ring to the list
411	export Jnew;
412	export eD;
413	export bM;
414	resList[size(resList)+1]=Rnew;
415	setring R;
416	kill Rnew;
417	}
418	setring S;
419	return(resList);
420	}
421	example
422	{ "EXAMPLE:";
423	echo = 2;
424	ring r=0,(x,y,z),dp;
425	ideal I=z2-x^3*y^2;
426	ideal C=z,xy;
427	list li=blowUp2(I,C);
428	size(li); // number of charts
429	def S1=li[1];
430	setring S1; // chart 1
431	basering;
432	Jnew;
433	eD;
434	bM;
435	def S2=li[2];
436	setring S2; // chart 2
437	basering;
438	Jnew;
439	eD;
440	bM;
441	}
442	///////////////////////////////////////////////////////////////////////////////
443
444	proc Center(ideal J,list #)
445	"USAGE: Center(J[,W][,E])
446	@* J,W = ideals
447	@* E = list
448	ASSUME: J = ideal containing W ( W = 0 if not specified)
449	@* E = list of smooth hypersurfaces (e.g. exceptional divisors)
450	COMPUTE: the center of the blow-up of J for the resolution algorithm
451	of [Bravo,Encinas,Villamayor]
452	RETURN: ideal, describing the center
453	EXAMPLE: example Center; shows an example
454	"
455	{
456	ideal W;
457	list E;
458	ideal abb=maxideal(1);
459	intvec v;
460	intvec bvec;
461	intvec w=-1;
462	matrix intE;
463	if(size(#)>0)
464	{
465	if(typeof(#[1])=="ideal")
466	{
467	W=#[1];
468	}
469	if(typeof(#[1])=="list")
470	{
471	E=#[1];
472	}
473	if(size(#)>1)
474	{
475	if((typeof(#[2])=="list") && (size(E)==0))
476	{
477	E=#[2];
478	}
479	if((typeof(#[2])=="ideal") && (size(W)==0))
480	{
481	W=#[2];
482	}
483	if(size(#)==3){bvec=#[3];}
484	}
485	}
486
487	list BO=W,J,bvec,E,abb,v,w,intE,intvec(0);
488	if(defined(invSat)){kill invSat;}
489	list invSat=ideal(0),intvec(0);
490	export(invSat);
491
492	list re=CenterBO(BO);
493	ideal cent=re[1];
494	return(cent);
495	}
496	example
497	{ "EXAMPLE:";
498	echo = 2;
499	ring R=0,(x,y),dp;
500	ideal J=x2-y3;
501	Center(J);
502	}
503	///////////////////////////////////////////////////////////////////////////////
504
505	proc blowUpBO(list BO, ideal C,int e)
506	"USAGE: blowUpBO (BO,C,e);
507	@* BO = basic object, a list: ideal W,
508	@* ideal J,
509	@* intvec b,
510	@* list Ex,
511	@* ideal ab,
512	@* intvec v,
513	@* intvec w,
514	@* matrix M
515	@* C = ideal
516	@* e = integer (0 usual blowing up, 1 deleting extra charts, 2 deleting
517	@* no charts )
518	ASSUME: R = basering, a polynomial ring, W an ideal of R,
519	@* J = ideal containing W,
520	@* C = ideal containing J
521	COMPUTE: the blowing up of BO[1] in C, the exeptional locus, the strict
522	transform of BO[2]
523	NOTE: blowUpBO may be applied to basic objects in the sense of
524	@* [Bravo, Encinas, Villamayor] in the following referred to as BO and
525	@* to presentations in the sense of [Bierstone, Milman] in the following
526	@* referred to as BM.
527	RETURN: a list l of length at most size(C),
528	l[i] is a ring containing an object BO resp. BM:
529	@* BO[1]=BM[1] an ideal, say Wi, defining the ambient space of the
530	i-th chart of the blowing up
531	@* BO[2]=BM[2] an ideal defining the strict transform
532	@* BO[3] intvec, the first integer b such that in the original object
533	(Delta^b(BO[2]))==1
534	the subsequent integers have the same property for Coeff-Objects
535	of BO[2] used when determining the center
536	@* BM[3] intvec, BM[3][i] is the assigned multiplicity of BM[2][i]
537	@* BO[4]=BM[4] the list of exceptional divisors
538	@* BO[5]=BM[5] an ideal defining the map K[W] ----> K[Wi]
539	@* BO[6]=BM[6] an intvec BO[6][j]=1 indicates that <BO[4][j],BO[2]>=1,
540	i.e. the strict transform does not meet the j-th exceptional
541	divisor
542	@* BO[7] intvec,
543	the index of the first blown-up object in the resolution process
544	leading to this object for which the value of b was BO[3]
545	the subsequent ones are the indices for the Coeff-Objects
546	of BO[2] used when determining the center
547	@* BM[7] intvec, BM[7][i] is the index at which the (2i-1)st entry
548	of the invariant first reached its current maximal value
549	@* BO[8]=BM[8] a matrix indicating that BO[4][i] meets BO[4][j] by
550	BO[8][i,j]=1 for i < j
551	@* BO[9] empty
552	@* BM[9] the invariant
553
554	EXAMPLE: example blowUpBO; shows an example
555	"
556	{
557	//---------------------------------------------------------------------------
558	// Initialization and sanity checks
559	//---------------------------------------------------------------------------
560	def R0=basering;
561	if(!defined(locaT)){int locaT;}
562	if(locaT){poly pp=@p;}
563	intvec v;
564	int shortC=defined(shortcut);
565	int invS=defined(invSat);
566	int eq,hy;
567	int extra,noDel,keepDiv;
568	if(e==1){extra=1;}
569	//---keeps all charts
570	if(e==2){noDel=1;}
571	//---this is only for curves and surfaces
572	//---keeps all charts with relevant informations on the exceptional divisors
573	if(e==3){keepDiv=1;}
574	if( typeof(attrib(BO[2],"isEqui"))=="int" )
575	{
576	eq=attrib(BO[2],"isEqui");
577	}
578	if( typeof(attrib(BO[2],"isHy"))=="int" )
579	{
580	hy=attrib(BO[2],"isHy");
581	}
582	string newvar;
583	int n=nvars(R0);
584	int i,j,l,m,x,jj,ll,haveCenters,co;
585	//---the center should neither be the whole space nor empty
586	if((size(C)==0)\|\|(deg(C[1])==0))
587	{
588	list result=R0;
589	return(result);
590	}
591	if(!defined(debugBlowUp))
592	{
593	int debugBlowUp=0;
594	}
595	//---------------------------------------------------------------------------
596	// Drop unnecessary variables
597	//---------------------------------------------------------------------------
598	//---step 1: substitution
599	if(!((keepDiv)\|\|(noDel)))
600	{
601	//!!! in case keepDiv and noDel:
602	//!!! maybe simplify the situation by an appropriate coordinate change
603	//!!! of this kind -- without dropping variables?
604	list L=elimpart(BO[1]);
605	if(size(L[2])!=0)
606	{
607	map psi=R0,L[5];
608	C=psi(C);
609	BO=psi(BO);
610	}
611
612	if(size(BO[1])==0)
613	{
614	ideal LL;
615	for(j=1;j<=size(BO[4]);j++)
616	{
617	LL=LL,BO[4][j];
618	}
619	LL=variables(LL);
620	L=elimpart(BO[2]);
621	if((size(L[2])!=0)&&(size(std(LL+L[2]))==size(L[2])+size(LL)))
622	{
623	map chi=R0,L[5];
624	C=chi(C);
625	BO=chi(BO);
626	}
627	}
628	//---step 2: dropping non-occurring variables
629	int s=size(C);
630	ideal K=BO[1],BO[2],C;
631	for(j=1;j<=size(BO[4]);j++)
632	{
633	K=K,BO[4][j];
634	}
635	list N=findvars(K);
636	if(size(N[1])<n)
637	{
638	newvar=string(N[1]);
639	v=N[4];
640	for(j=1;j<=size(v);j++){BO[5]=subst(BO[5],var(v[j]),0);}
641	execute("ring R1=("+charstr(R0)+"),("+newvar+"),dp;");
642	list BO=imap(R0,BO);
643	ideal C=imap(R0,C);
644	n=nvars(R1);
645	}
646	else
647	{
648	def R1=basering;
649	}
650	}
651	else
652	{
653	int s=size(C);
654	def R1=basering;
655	}
656	if(debugBlowUp)
657	{
658	"---> In BlowUp: After dropping unnecessary variables";
659	"BO:";
660	BO;
661	"C:";
662	C;
663	}
664	//---------------------------------------------------------------------------
665	// Do the actual blow-up
666	//---------------------------------------------------------------------------
667	//--- control the names of the variables
668	execute("ring R=("+charstr(R0)+"),(x(1..n)),dp;");
669	list BO=fetch(R1,BO);
670	ideal C=fetch(R1,C);
671	list Cmstd=mstd(C);
672	C=Cmstd[2];
673	if(size(Cmstd[1])<=size(Cmstd[2]))
674	{
675	C=Cmstd[1];
676	}
677	else
678	{
679	C=interred(C);
680	}
681	list result;
682	//--- the blow-up process
683	ideal W =BO[1];
684	ideal J =BO[2];
685	intvec bvec =BO[3];
686	list Ex =BO[4];
687	ideal abb=BO[5];
688	intvec wvec=BO[7];
689	ideal laM=maxideal(1);
690	if((typeof(BO[9])=="intmat")\|\|(typeof(BO[9])=="intvec"))
691	{
692	def @invmat=BO[9];
693	}
694	if(size(BO)>9)
695	{
696	//--- check whether a previous center had been split into connected components
697	if(size(BO[10])>0)
698	{
699	list knownCenters=BO[10];
700	haveCenters=1;
701	}
702	}
703
704	matrix intE=BO[8];
705	Ex[size(Ex)+1]=var(1);
706	//to have the list depending on R in case BO[4] is empty
707
708	execute("ring S=("+charstr(R)+"),("+varstr(R)+",y(0..s-1)),dp;");
709	list resu;
710	list B;
711	execute("ring T=("+charstr(R)+"),("+varstr(R)+",y(0..s-1),t),dp;");
712	ideal C=imap(R,C);
713	ideal W=imap(R,W);
714	execute("map phi=S,"+varstr(R)+",t*C;")
715	setring S;
716
717	//--- the ideal describing the blow-up map
718	ideal abb=imap(R,abb);
719	ideal laM0=imap(R,laM);
720	//--- the ideal of the blowing up of the ambient space
721	ideal W=preimage(T,phi,W);
722	//--- the ideal of the exceptional locus
723	ideal E=imap(R,C);
724
725	list E1=imap(R,Ex);
726	E1[size(E1)]=E;
727	ideal J=imap(R,J)+W;
728
729	if(haveCenters){list kN=imap(R,knownCenters);}
730
731	//--- the strict transform of the exceptional divisors
732	for(j=1;j<size(E1);j++){E1[j]=sat(E1[j]+W,E)[1];}
733	//--- the intersection matrix of the exceptional divisors
734	matrix intEold=imap(R,intE);
735	matrix intE[size(E1)][size(E1)];
736	ideal U,Jsub,sLstd;
737
738	for(j=1;j<size(E1);j++)
739	{
740	for(l=j+1;l<=size(E1);l++)
741	{
742	if(deg(E1[j][1])==0)
743	{
744	if(l<size(E1)){intE[j,l]=intEold[j,l];}
745	else {intE[j,l]=0;}
746	}
747	else
748	{
749	if(deg(E1[l][1])==0)
750	{
751	if(l<size(E1)){intE[j,l]=intEold[j,l];}
752	}
753	else
754	{
755	U= std(E1[l]+E1[j]);
756	if(dim(U)>0){intE[j,l]=1;}
757	else {intE[j,l]=0;}
758	}
759	}
760	}
761	}
762	if(debugBlowUp)
763	{
764	"----> In BlowUp: After Blowing-up, before Clean-Up";
765	"W:";
766	W;
767	"J:";
768	J;
769	}
770	//----------------------------------------------------------------------------
771	// generating and cleaning up the different charts
772	//----------------------------------------------------------------------------
773	list M;
774	map psi;
775	list E2;
776	ideal K,JJ,laM,LL,MM;
777	n=nvars(S);
778	list N;
779	list Bstd;
780	intvec delCharts,extraCharts;
781	delCharts[s]=0;
782	extraCharts[s]=0;
783	ideal MA=y(0..s-1);
784	list ZRes,ZlaM,ZsLstd;
785	for(i=0;i<=s-1;i++)
786	{
787	if(haveCenters)
788	{
789	B[10]=kN;
790	for(j=1;j<=size(kN);j++)
791	{
792	B[10][j][1]=subst(B[10][j][1],y(i),1);
793	}
794	}
795	B[8]=intE;
796	B[1]=std(subst(W,y(i),1));
797	if(deg(B[1][1])==0)
798	{
799	//--- subsets of the empty set are not really interesting!
800	delCharts[i+1]=1;
801	ZRes[i+1]=B;
802	ZlaM[i+1]=laM;
803	i++;
804	continue;
805	}
806	Jsub=subst(J,y(i),1);
807	attrib(Jsub,"isEqui",eq);
808	attrib(Jsub,"isHy",hy);
809	B[2]=Jsub;
810	B[3]=bvec;
811	for(j=1;j<size(E1);j++){E2[j]=subst(E1[j],y(i),1);}
812	E2[size(E1)]=E+B[1];
813	B[4]=E2;
814	M=elimpart(B[1]);
815	B[5]=abb;
816	laM=laM0;
817	psi=S,maxideal(1);
818	if(size(M[2])!=0)
819	{
820	psi=S,M[5];
821	B=psi(B);
822	laM=psi(laM);
823	}
824	Jsub=B[2];
825	B[2]=sat(Jsub,std(psi(E)))[1];
826	if(!defined(MAtmp)){ideal MAtmp=MA;}
827	MAtmp[i+1]=0;
828	JJ=std(B[2]+MAtmp);
829	if(deg(JJ[1])==0)
830	{
831	delCharts[i+1]=1;
832	//--- the i-th chart will be marked for deleting because all informations
833	//--- are already contained in the union of the remaining charts
834	}
835	else
836	{
837	if((eq)&&(dim(JJ)<dim(std(B[2]))))
838	{
839	extraCharts[i+1]=1;
840	//--- compute the singular locus
841	if((B[3][1]<=1)&&(hy))
842	//--- B[2] is a smooth hypersurface
843	{
844	sLstd=ideal(1);
845	}
846	else
847	{
848	Bstd=mstd(B[2]);
849	if(n-dim(Bstd[1])>4)
850	{
851	//--- in this case the singular locus is too complicated
852	sLstd=ideal(0);
853	}
854	JJ=Bstd[2];
855	attrib(JJ,"isEqui",eq);
856	B[2]=JJ;
857	sLstd=slocusE(B[2]);
858	}
859	m=0;
860	if(deg(std(sLstd+MAtmp)[1])==0)
861	{
862	//--- the singular locus of B[2] is in the union of the remaining charts
863	m=1;
864	for(l=1;l<=size(B[4]);l++)
865	{
866	if(deg(std(B[2]+B[4][l]+MAtmp)[1])!=0)
867	{
868	//--- the exceptional divisor meets B[2] at the locus of MAtmp
869	//--- we continue only if the option extra=1 and we have transversal
870	//--- intersection
871	m=0;
872	break;
873	}
874	}
875	}
876	if(m)
877	{
878	//--- the i-th chart will be marked for deleting because all informations
879	//--- are already contained in the union of the remaining charts
880	delCharts[i+1]=1;
881	}
882	}
883	}
884	if(delCharts[i+1]==0)
885	{
886	MAtmp[i+1]=MA[i+1];
887	ZsLstd[i+1]=sLstd;
888	}
889	ZRes[i+1]=B;
890	ZlaM[i+1]=laM;
891	}
892	//---------------------------------------------------------------------------
893	// extra = ignore uninteresting charts even if there is a normal
894	// crosssing intersection in it
895	//---------------------------------------------------------------------------
896	if(extra)
897	{
898	for(i=0;i<=s-1;i++)
899	{
900	if((delCharts[i+1]==0)&&(extraCharts[i+1]))
901	{
902	MAtmp[i+1]=0;
903	B=ZRes[i+1];
904	sLstd=ZsLstd[i+1];
905	m=0;
906	if(deg(std(sLstd+MAtmp)[1])==0)
907	{
908	//--- the singular locus of B[2] is in the union of the remaining charts
909	m=1;
910	for(l=1;l<=size(B[4]);l++)
911	{
912	if(deg(std(B[2]+B[4][l]+MAtmp)[1])!=0)
913	{
914	//--- the exceptional divisor meets B[2] at the locus of MAtmp
915	//--- we continue only if the option extra=1 and we have transversal
916	//--- intersection
917	m=2;
918	if(!transversalTB(B[2],list(B[4][l]),MAtmp))
919	{
920	m=0;break;
921	}
922	}
923	}
924	}
925	if(m)
926	{
927	if(m==1)
928	{
929	//--- the i-th chart will be marked for deleting because all informations
930	//--- are already contained in the union of the remaining charts
931	delCharts[i+1]=1;
932	}
933	else
934	{
935	//--- the option extra=1 and we have transversal intersection
936	//--- we delete the chart in case of normal crossings
937	if(normalCrossB(B[2],B[4],MAtmp))
938	{
939	//--- in case of the option extra
940	//--- the i-th chart will be marked for deleting because all informations
941	//--- are already contained in the union of the remaining charts
942
943	delCharts[i+1]=1;
944	}
945	}
946	}
947	if(delCharts[i+1]==0)
948	{
949	MAtmp[i+1]=MA[i+1];
950	}
951	}
952	}
953	for(i=0;i<=s-1;i++)
954	{
955	if(!delCharts[i+1]){break;}
956	}
957	if(i==s){delCharts[s]=0;}
958	}
959	for(i=0;i<=s-1;i++)
960	{
961	B=ZRes[i+1];
962	laM=ZlaM[i+1];
963	if(noDel){delCharts[i+1]=0;}
964	//---keeps chart if the exceptional divisor is not in any other chart
965	if((delCharts[i+1])&&(keepDiv))
966	{
967	for(j=1;j<=size(B[4]);j++)
968	{
969	if(deg(std(B[4][j])[1])>0)
970	{
971	x=0;
972	for(l=0;l<=s-1;l++)
973	{
974	if((l!=i)&&(!delCharts[l+1])&&(deg(std(ZRes[l+1][4][j])[1])>0))
975	{
976	x=1;
977	break;
978	}
979	}
980	if(!x)
981	{
982	delCharts[i+1]=0;
983	//!!!evtl. diese Karten markieren und nicht weiter aufblasen???
984	break;
985	}
986	}
987	}
988	}
989	//---keeps charts if the intersection of 2 divisors is not in any other chart
990	if((delCharts[i+1])&&(keepDiv))
991	{
992	for(j=1;j<=size(B[4])-1;j++)
993	{
994	for(l=j+1;l<=size(B[4]);l++)
995	{
996	if(deg(std(B[4][j]+B[4][l])[1])>0)
997	{
998	x=0;
999	for(ll=0;ll<=s-1;ll++)
1000	{
1001	if((ll!=i)&&(!delCharts[ll+1]))
1002	{
1003	if(deg(std(ZRes[ll+1][4][j]+ZRes[ll+1][4][l])[1])>0)
1004	{
1005	x=1;
1006	break;
1007	}
1008	}
1009	}
1010	if(!x)
1011	{
1012	delCharts[i+1]=0;
1013	break;
1014	}
1015	}
1016	}
1017	if(!delCharts[i+1]){break;}
1018	}
1019	}
1020	//---keeps charts if the intersection of 3 divisors is not in any other chart
1021	if((delCharts[i+1])&&(keepDiv))
1022	{
1023	for(j=1;j<=size(B[4])-2;j++)
1024	{
1025	for(l=j+1;l<=size(B[4])-1;l++)
1026	{
1027	for(ll=l+1;ll<=size(B[4]);ll++)
1028	{
1029	if(deg(std(B[4][j]+B[4][l]+B[4][ll])[1])>0)
1030	{
1031	x=0;
1032	for(jj=0;jj<=s-1;jj++)
1033	{
1034	if((jj!=i)&&(!delCharts[jj+1]))
1035	{
1036	if(deg(std(ZRes[jj+1][4][j]
1037	+ZRes[jj+1][4][l]+ZRes[jj+1][4][ll])[1])>0)
1038	{
1039	x=1;
1040	break;
1041	}
1042	}
1043	}
1044	if(!x)
1045	{
1046	delCharts[i+1]=0;
1047	break;
1048	}
1049	}
1050	}
1051	if(!delCharts[i+1]){break;}
1052	}
1053	if(!delCharts[i+1]){break;}
1054	}
1055	}
1056	if(delCharts[i+1]==0)
1057	{
1058	//--- try to decrease the number of variables by substitution
1059	if((!keepDiv)&&(!noDel))
1060	{
1061	list WW=elimpart(B[1]);
1062	map phiW=basering,WW[5];
1063	B=phiW(B);
1064	laM=phiW(laM);
1065	kill WW;
1066	kill phiW;
1067	if(size(B[1])==0)
1068	{
1069	LL=0;
1070	for(j=1;j<=size(B[4]);j++)
1071	{
1072	MM=std(B[4][j]);
1073	if(deg(MM[1])>0){LL=LL,MM;}
1074	}
1075	LL=variables(LL);
1076	M=elimpart(B[2]);
1077	if((size(M[2])!=0)&&(size(std(LL+M[2]))==size(M[2])+size(LL)))
1078	{
1079	psi=S,M[5];
1080	B=psi(B);
1081	laM=psi(laM);
1082	}
1083	}
1084	}
1085	//---- interreduce B[1],B[2] and all B[4][j]
1086	B[1]=interred(B[1]);
1087	B[2]=interred(B[2]);
1088	E2=B[4];
1089	for(j=1;j<=size(E2);j++){E2[j]=interred(E2[j]);}
1090	B[4]=E2;
1091	v=0;v[size(E2)]=0;
1092	//--- mark those j for which B[4] does not meet B[2]
1093	for(j=1;j<=size(E2);j++)
1094	{
1095	K=E2[j],B[2];
1096	K=std(K);
1097	if(deg(K[1])==0)
1098	{
1099	v[j]=1;
1100	}
1101	}
1102	B[6]=v;
1103	B[7]=wvec;
1104	//--- throw away variables which do not occur
1105	K=B[1],B[2],B[5]; //Aenderung!!!
1106	for(j=1;j<=size(B[4]);j++){K=K,B[4][j];}
1107	N=findvars(K);
1108	if(size(N[1])<n)
1109	{
1110	newvar=string(N[1]);
1111	v=N[4];
1112	for(j=1;j<=size(v);j++)
1113	{
1114	B[5]=subst(B[5],var(v[j]),0);
1115	laM=subst(laM,var(v[j]),0);
1116	}
1117	execute("ring R2=("+charstr(S)+"),("+newvar+"),dp;");
1118	list BO=imap(S,B);
1119	ideal laM=imap(S,laM);
1120	}
1121	else
1122	{
1123	def R2=basering;
1124	list BO=B;
1125	}
1126	ideal JJ=BO[2];
1127	attrib(JJ,"isEqui",eq);
1128	attrib(JJ,"isHy",hy);
1129	BO[2]=JJ;
1130	//--- strict transforms of the known centers
1131	if(haveCenters)
1132	{
1133	ideal tt;
1134	list tList;
1135	for(j=1;j<=size(BO[10]);j++)
1136	{
1137	tt=std(BO[10][j][1]);
1138	if(deg(tt[1])>0)
1139	{
1140	tt=sat(tt,BO[4][size(BO[4])])[1];
1141	}
1142	if((deg(tt[1])>0)&&(deg(std(tt+BO[2]+BO[1])[1])>0))
1143	{
1144	tList[size(tList)+1]=
1145	list(tt,BO[10][j][2],BO[10][j][3],BO[10][j][4]);
1146	}
1147	}
1148	BO[10]=tList;
1149	kill tList;
1150	}
1151	//--- marking variables which do not occur in BO[1] and BO[2]
1152	//--- and occur in exactly one BO[4][j], which is the hyperplane given by
1153	//--- this variable
1154	//!!!! not necessarily in exactly one BO[4][j]
1155	list N=findvars(BO[1]+BO[2]);
1156	if(size(N[1])<nvars(basering))
1157	{
1158	v=N[4];
1159	if(defined(H)){kill H;}
1160	if(defined(EE)){kill EE;}
1161	if(defined(vv)){kill vv;}
1162	list EE;
1163	intvec vv;
1164	ideal H=maxideal(1);
1165	for(j=1;j<=size(v);j++)
1166	{
1167	H[v[j]]=0;
1168	}
1169	H=std(H);
1170	for(l=1;l<=size(BO[4]);l++)
1171	{
1172	if(BO[6][l]==1){l++;continue;}
1173	if(size(ideal(matrix(reduce(BO[4][l],H))-BO[4][l]))==0)
1174	{
1175	//!!! need further cleanup:
1176	//!!! this part of the code is no longer used since it did not glue properly
1177	// BO[6][l]=2;
1178	EE[size(EE)+1]=BO[4][l];
1179	vv[size(vv)+1]=l;
1180	}
1181	}
1182	if((size(vv)>dim(std(BO[2])))&&(deg(BO[2][1])>0))
1183	{
1184	list BOtemp3=BO;
1185	BOtemp3[4]=EE;
1186	intvec ivtemp3;
1187	ivtemp3[size(BOtemp3[4])]=0;
1188	BOtemp3[6]=ivtemp3;
1189	BOtemp3[7][1]=-1;
1190	list iEtemp3=inters_E(BOtemp3);
1191	if(iEtemp3[2]>=dim(std(BOtemp3[2])))
1192	{
1193	for(l=2;l<=size(vv);l++)
1194	{
1195	BO[6][vv[l]]=0;
1196	}
1197	}
1198	kill BOtemp3,ivtemp3,iEtemp3;
1199	}
1200	}
1201	list thisChart=ideal(0),i;
1202	export thisChart;
1203
1204	//----------------------------------------------------------------------------
1205	// export the basic object and append the ring to the list of rings
1206	//----------------------------------------------------------------------------
1207	if(debugBlowUp)
1208	{
1209	"----> In BlowUp: Adding a single chart";
1210	"BO:";
1211	BO;
1212	}
1213	if(locaT)
1214	{
1215	map locaPhi=R0,laM;
1216	poly @p=locaPhi(pp);
1217	export(@p);
1218	}
1219	ideal lastMap=laM;
1220	export lastMap;
1221	if(invS){list invSat=imap(R0,invSat);export invSat;}
1222	if(defined(@invmat)){BO[9]=@invmat;}
1223	if(shortC){list shortcut=imap(R0,shortcut);export(shortcut);}
1224	export BO;
1225	result[size(result)+1]=R2;
1226	setring S;
1227	kill R2;
1228	}
1229	}
1230	setring R0;
1231	return(result);
1232	}
1233	example
1234	{"EXAMPLE:";
1235	echo = 2;
1236	ring R=0,(x,y),dp;
1237
1238	ideal W;
1239	ideal J=x2-y3;
1240	intvec b=1;
1241	list E;
1242	ideal abb=maxideal(1);
1243	intvec v;
1244	intvec w=-1;
1245	matrix M;
1246	intvec ma;
1247	list BO=W,J,b,E,abb,v,w,M,ma;
1248
1249	ideal C=CenterBO(BO)[1];
1250
1251	list blow=blowUpBO(BO,C,0);
1252	def Q=blow[1];
1253	setring Q;
1254	BO;
1255	}
1256	//////////////////////////////////////////////////////////////////////////////
1257	static proc slocusE(ideal i)
1258	"Internal procedure - no help and no example available
1259	"
1260	{
1261	//--- do slocus in equidimensional case directly -- speed up
1262	if(size(i)==0){return(ideal(1));}
1263	if( typeof(attrib(i,"isEqui"))=="int" )
1264	{
1265	if(attrib(i,"isEqui")==1)
1266	{
1267	ideal j=std(i);
1268	if(deg(j[1])==0){return(j);}
1269	int cod = nvars(basering)-dim(j);
1270	i = i+minor(jacob(i),cod);
1271	return(i);
1272	}
1273	}
1274	return(slocus(i));
1275	}
1276	///////////////////////////////////////////////////////////////////////////////
1277	static proc inters_E(list BO)
1278	"USAGE: inters_E(BO);
1279	@* BO = basic object, a list: ideal W,
1280	@* ideal J,
1281	@* intvec b,
1282	@* list Ex,
1283	@* ideal ab,
1284	@* intvec v,
1285	@* intvec w,
1286	@* matrix M
1287	ASSUME: R = basering, a polynomial ring, W an ideal of R,
1288	@* J = ideal containing W,
1289	@* BO in the setting of case 2 of [Bravo,Encinas,Villamayor]
1290	@* BO[4]=E, BO[4][1..count]=E^-
1291	@* BO[7][1]=count
1292	COMPUTE: (W,(P,1),E^+) in the notation of [Bravo,Encinas,Villamayor]
1293	RETURN: a list l ,
1294	l[1]: P = product of ideals I(H_i1)+..+I(H_in) over all
1295	n-tuples of indices i1..in from 1..count
1296	l[2]: n = maximal number of H_i from E^- meeting J simultaneously
1297	l[3]: maximal number of H_i from E meeting J simultaneously
1298	EXAMPLE: internal procedure - no example available
1299	"
1300	{
1301	//---------------------------------------------------------------------
1302	// Initialization
1303	//---------------------------------------------------------------------
1304	int kk,jj,ii,updated,n,count2,kkdiff;
1305	def rb=basering;
1306	def W=BO[1];
1307	ideal J=BO[1],BO[2];
1308	int nonnormal;
1309	int maxkk=dim(std(J));
1310	int dimJ=maxkk;
1311	ideal test2;
1312	list merklist1,merklist2;
1313	if(size(BO[4])==0)
1314	{
1315	list retlist=BO[2],n;
1316	return(retlist);
1317	}
1318	def E=BO[4];
1319	intvec stoplist=BO[6];
1320	//--- fill in all known information about exceptional divisors not meeting
1321	//--- current chart
1322	for(ii=1;ii<=size(E);ii++)
1323	{
1324	if(deg(std(E[ii])[1])==0)
1325	{
1326	stoplist[ii]=1;
1327	}
1328	}
1329
1330	int count=BO[7][1];
1331	if(!defined(debug_Inters_E))
1332	{
1333	int debug_Inters_E=0;
1334	}
1335	//---------------------------------------------------------------------
1336	// we only want to look at E^-, not at all of E
1337	//---------------------------------------------------------------------
1338	if (count>-1)
1339	{
1340	if (count>0)
1341	{
1342	list E_new=E[1..count];
1343	count2=size(E);
1344	}
1345	else
1346	{
1347	list E_new;
1348	count2=size(E);
1349	}
1350	}
1351	else
1352	{
1353	list E_new=E;
1354	count=size(E_new);
1355	count2=count;
1356	}
1357	//---------------------------------------------------------------------
1358	// combinatorics is expensive in an interpreted language,
1359	// we leave it to the kernel by translating it into monomial
1360	// ideals in a new ring with variables t(i)
1361	//---------------------------------------------------------------------
1362	string rstr="ring rcomb=0,(t(1.." + string(size(E)) + ")),dp;";
1363	execute(rstr);
1364	ideal potid,potid2;
1365	list monlist,comblist,merkmon;
1366	for(kk=1;kk<=count;kk++)
1367	{
1368	if(stoplist[kk]==0)
1369	{
1370	//**************************************************************************/
1371	// it does not make sense to intersect twice by the same E_i
1372	// ===> reduce by t(i)^2
1373	//**************************************************************************/
1374	potid=potid,t(kk)^2;
1375	}
1376	else
1377	{
1378	//**************************************************************************/
1379	// it does not make sense to consider E_i with v[i]==1 ===> reduce by t(i)
1380	//**************************************************************************/
1381	potid=potid,t(kk);
1382	//**************************************************************************/
1383	// if stoplist[kk]==2 then J and all E_i automatically intersect E_kk
1384	// hence we need not test it, but we have to lower maxkk by one
1385	//**************************************************************************/
1386	if(stoplist[kk]==2)
1387	{
1388	maxkk=maxkk-1;
1389	kkdiff++; // count these for dimension check later on
1390
1391	}
1392	}
1393	}
1394
1395	potid2=std(potid);
1396	if(count2>count)
1397	{
1398	potid=potid,t((count+1)..count2);
1399	for(kk=max(1,count);kk<=count2;kk++)
1400	{
1401	potid2=potid2,t(kk)^2;
1402	}
1403	}
1404	potid=std(potid);
1405	potid2=std(potid2);
1406	for(kk=1;(((kk<=count)\|\|(kk<=maxkk+1))&&(kk<=count2));kk++)
1407	{
1408	//-------------------------------------------------------------------------
1409	// monlist[kk]=lists of kk entries of E_new, not containing an E_i twice,
1410	// not containing an E_i where v[i]==1
1411	//-------------------------------------------------------------------------
1412	monlist[kk]=redMax(kk,potid);
1413	//*************************************************************************/
1414	// in the case of n<=maxkk we also need to know whether n would still be
1415	// below this bound if we considered all of E instead of E_new
1416	// ===> merkmon contains previously ignored tuples E_i1,..,E_im
1417	//*************************************************************************/
1418	if(kk<=maxkk+1)
1419	{
1420	merkmon[kk]=redMax(kk,potid2);
1421	merkmon[kk]=simplify(reduce(merkmon[kk],std(monlist[kk])),2);
1422	}
1423	}
1424	if(debug_Inters_E)
1425	{
1426	"----> In Inters_E: the tuples";
1427	"tuples of E^-:";
1428	monlist;
1429	"the remaining tuples:";
1430	merkmon;
1431	}
1432	//-------------------------------------------------------------------------
1433	// check whether there is a kk-tuple of E_i intersecting J,
1434	// kk running from 1 to count
1435	//-------------------------------------------------------------------------
1436	for(kk=1;kk<=count;kk++)
1437	{
1438	if(size(monlist[kk])==0) break;
1439	kill comblist;
1440	list comblist;
1441	//--- transscribe the tuples from monomial notation to intvec notation
1442	for(jj=1;jj<=ncols(monlist[kk]);jj++)
1443	{
1444	comblist[jj]=leadexp(monlist[kk][jj]);
1445	}
1446	setring rb;
1447	updated=0;
1448	//------------------------------------------------------------------------
1449	// Do the intersections
1450	//------------------------------------------------------------------------
1451	for(jj=1;jj<=size(comblist);jj++)
1452	{
1453	//--- jj-th tuple from list of tuples of kk E_i
1454	test2=J;
1455	for(ii=1;ii<=count;ii++)
1456	{
1457	if(comblist[jj][ii]==1)
1458	{
1459	test2=test2,E_new[ii];
1460	}
1461	}
1462	test2=std(test2);
1463	//--- check whether this intersection is non-empty and store it accordingly
1464	if(deg(test2[1])!=0)
1465	{
1466	//--- it is non-empty
1467	if(updated!=0)
1468	{
1469	merklist1[size(merklist1)+1]=comblist[jj];
1470	}
1471	else
1472	{
1473	kill merklist1;
1474	list merklist1;
1475	merklist1[1]=comblist[jj];
1476	updated=1;
1477	n=kk;
1478	}
1479	if(dim(test2)!=maxkk-kk+kkdiff)
1480	{
1481	nonnormal=1;
1482	}
1483	}
1484	else
1485	{
1486	//--- it is empty
1487	merklist2[size(merklist2)+1]=jj;
1488	}
1489	}
1490	setring rcomb;
1491	ideal redid;
1492	//---------------------------------------------------------------------
1493	// update monlist and merkmon by the knowledge what intersections are
1494	// empty in the kk-th step
1495	//---------------------------------------------------------------------
1496	for(jj=1;jj<=size(merklist2);jj++)
1497	{
1498	redid=redid,monlist[kk][merklist2[jj]];
1499	}
1500	for(jj=kk+1;jj<=count;jj++)
1501	{
1502	monlist[jj]=simplify(reduce(monlist[jj],std(redid)),2);
1503	if(jj<=maxkk+1)
1504	{
1505	merkmon[jj]=simplify(reduce(merkmon[jj],std(redid)),2);
1506	}
1507	}
1508	kill redid;
1509	kill merklist2;
1510	list merklist2;
1511	}
1512	if(debug_Inters_E)
1513	{
1514	"----> In Inters_E: intersections found:";
1515	merklist1;
1516	}
1517	//---------------------------------------------------------------------
1518	// form the union of the intersections of the appropriate E_i
1519	//---------------------------------------------------------------------
1520	setring rb;
1521	ideal center,dummy;
1522	list centlist;
1523	for(kk=1;kk<=size(merklist1);kk++)
1524	{
1525	for(jj=1;jj<=size(merklist1[kk]);jj++)
1526	{
1527	if(merklist1[kk][jj]==1)
1528	{
1529	dummy=dummy,E_new[jj];
1530	}
1531	}
1532	if(size(center)==0)
1533	{
1534	center=dummy;
1535	centlist[1]=dummy;
1536	}
1537	else
1538	{
1539	center=intersect(center,dummy);
1540	centlist[size(centlist)+1]=dummy;
1541	}
1542	dummy=0;
1543	}
1544	if(debug_Inters_E)
1545	{
1546	"----> In Inters_E: intersection of E_i";
1547	"maximal number of E_i encountered in:";
1548	center;
1549	"the components of this locus:";
1550	centlist;
1551	"maximal number of E_i from E^- intersecting simultaneously:",n;
1552	if(nonnormal)
1553	{
1554	"flag nonnormal is set";
1555	}
1556	}
1557	list retlist=center,n;
1558	//-------------------------------------------------------------------------
1559	// If n<=maxkk, then test if this is the case for all of E not just E_new
1560	// using the pairs indicated by merkmon
1561	//-------------------------------------------------------------------------
1562	int ntotal=n;
1563	if((n<=maxkk)&&(n<count2)&&(!nonnormal))
1564	{
1565	//--- check kk-tuples
1566	for(kk=n+1;kk<=maxkk+1;kk++)
1567	{
1568	setring rcomb;
1569	//--- check if there are combinations to be checked
1570	if(kk>size(merkmon))
1571	{
1572	setring rb;
1573	break;
1574	}
1575	if(size(merkmon[kk])!=0)
1576	{
1577	kill comblist;
1578	list comblist;
1579	//--- transscribe tuples from monomial notation to intvec notation
1580	for(jj=1;jj<=size(merkmon[kk]);jj++)
1581	{
1582	comblist[jj]=leadexp(merkmon[kk][jj]);
1583	}
1584	setring rb;
1585	//--- check jj-th tuple from the list of kk-tuples
1586	for(jj=1;jj<=size(comblist);jj++)
1587	{
1588	test2=J;
1589	for(ii=1;ii<=nvars(rcomb);ii++)
1590	{
1591	if(comblist[jj][ii]==1)
1592	{
1593	test2=test2,E[ii];
1594	}
1595	}
1596	test2=std(test2);
1597	//--- as soon as we found one we can proceed to the subsequent kk
1598	if(deg(test2[1])!=0)
1599	{
1600	ntotal=kk;
1601	if(dim(test2)-kkdiff!=maxkk-kk)
1602	{
1603	nonnormal=2;
1604	break;
1605	}
1606	}
1607	}
1608	//--- if we already know that too many E_i intersect simultaneously,
1609	//--- we need not proceed any further
1610	if(nonnormal)
1611	{
1612	break;
1613	}
1614	}
1615	else
1616	{
1617	setring rb;
1618	break;
1619	}
1620	}
1621	}
1622	//-------------------------------------------------------------------------
1623	// update the result accordingly and return it
1624	//-------------------------------------------------------------------------
1625	if(maxkk<dimJ)
1626	{
1627	n=n+dimJ-maxkk;
1628	ntotal=ntotal+dimJ-maxkk;
1629	}
1630	retlist[2]=n;
1631	retlist[3]=ntotal;
1632	if(n<=dimJ)
1633	{
1634	retlist[4]=centlist;
1635	retlist[5]=merklist1;
1636	if(nonnormal)
1637	{
1638	retlist[6]=nonnormal;
1639	}
1640	}
1641	return(retlist);
1642	}
1643	///////////////////////////////////////////////////////////////////////////
1644
1645	proc Delta(list BO)
1646	"USAGE: Delta (BO);
1647	@* BO = basic object, a list: ideal W,
1648	@* ideal J,
1649	@* intvec b,
1650	@* list Ex,
1651	@* ideal ab,
1652	@* intvec v,
1653	@* intvec w,
1654	@* matrix M
1655	ASSUME: R = basering, a polynomial ring, W an ideal of R,
1656	@* J = ideal containing W
1657	COMPUTE: Delta-operator applied to J in the notation of
1658	[Bravo,Encinas,Villamayor]
1659	RETURN: ideal
1660	EXAMPLE: example Delta; shows an example
1661	"
1662	{
1663	//---------------------------------------------------------------------------
1664	// Initialization and sanity checks
1665	//---------------------------------------------------------------------------
1666	ideal W=BO[1];
1667	ideal J=BO[2];
1668	ideal C=simplify(reduce(J,std(W)),2);
1669	list LC;
1670	int n=nvars(basering);
1671	//---------------------------------------------------------------------------
1672	// Simple case: W is the empty set
1673	//---------------------------------------------------------------------------
1674	if(size(W)==0)
1675	{
1676	C=C,jacob(J);
1677	C=std(C);
1678	return(C);
1679	}
1680	//---------------------------------------------------------------------------
1681	// General case: W is non-empty
1682	// Step 1: Find a minor of the Jacobian of W which is not identically zero
1683	// and look at the complement of the zero-set given by this minor;
1684	// this leads to the system of local parameters
1685	// Step 2: Form the derivatives w.r.t. this system of parameters
1686	//---------------------------------------------------------------------------
1687	//--- Step 1
1688	list re=findMinor(W);
1689	list L;
1690	int ii,i,j,l,k;
1691	J=C;
1692	ideal D=ideal(1);
1693	intvec v,w;
1694	ideal V;
1695	poly m;
1696
1697	for(ii=1;ii<=size(re);ii++)
1698	{
1699	C=J;
1700	L=re[ii];
1701	matrix A=L[1]; //(1/m)*A is the inverse matrix of the Jacobian of W
1702	//corresponding to m
1703	m=L[2]; //a k- minor of jacob(W), not identically zero
1704	//k=n-dim(W)
1705	V=L[3]; //the elements of W corresponding to m
1706	v=L[4]; //the indices of variables corresponding to m
1707	w=L[5]; //the indices of variables not corresponding to m
1708
1709	//--- Step 2
1710	//--- first some initialization depending on results of step 1
1711	k=size(V);
1712	matrix dg[1][k];
1713	matrix df[k][1];
1714	//--- derivatives of the generators of J w.r.t. system of parameters
1715	for(i=1;i<=size(J);i++)
1716	{
1717	for(j=1;j<=n-k;j++)
1718	{
1719	for(l=1;l<=k;l++)
1720	{
1721	dg[1,l]=diff(V[l],var(w[j]));
1722	df[l,1]=diff(J[i],var(v[l]));
1723	}
1724	C=C,mdiff(J[i],var(w[j]))-dgA*df;
1725	}
1726	}
1727	//--- everything should live in W, not just in the intersection of
1728	//--- D(m) with W
1729	C=C+W;
1730	C=sat(C,m)[1];
1731	//--- intersect ideal with previously computed ones to make sure that no
1732	//--- components are lost
1733	D=intersect(D,C);
1734	kill dg,df,A;
1735	}
1736	//--- return minimal set of generators of the result
1737	list li=mstd(D);
1738	D=li[2];
1739	if(size(li[1])<=size(D)){D=li[1];}
1740	return(D);
1741	}
1742	example
1743	{ "EXAMPLE:";
1744	echo = 2;
1745	ring R=0,(x,y,z),dp;
1746
1747	ideal W=z^2-x;
1748	ideal J=x*y^2+x^3;
1749	intvec b=1;
1750	list E;
1751	ideal abb=maxideal(1);
1752	intvec v;
1753	intvec w=-1;
1754	matrix M;
1755
1756	list BO=W,J,b,E,abb,v,w,M;
1757
1758	Delta(BO);
1759	}
1760
1761	//////////////////////////////////////////////////////////////////////////////
1762	static proc redMax(int k,ideal J)
1763	"Internal procedure - no help and no example available
1764	"
1765	{
1766	//--- reduce maxideal(k) by J, more efficient approach
1767	int i;
1768	ideal K=simplify(reduce(maxideal(1),J),2);
1769	for(i=2;i<=k;i++)
1770	{
1771	K=simplify(reduce(K*maxideal(1),J),2);
1772	}
1773	return(K);
1774	}
1775
1776	//////////////////////////////////////////////////////////////////////////////
1777	static proc findMinor(ideal W)
1778	"Internal procedure - no help and no example available
1779	"
1780	{
1781	//---------------------------------------------------------------------------
1782	// Initialization and sanity checks
1783	//---------------------------------------------------------------------------
1784	list L;
1785	intvec v,w;
1786	ideal Wstd=std(W);
1787	int n=nvars(basering); // total number of columns of Jacobian
1788	int k=n-dim(Wstd); // size of minors of Jacobian
1789	int a=size(W); // total number of rows of Jacobian
1790	matrix A[k][k];
1791	list LW=indexSet(a,k); // set of tuples of k rows
1792	list LV=indexSet(n,k); // set of tuples of k columns
1793	ideal IW,IV;
1794	int i,j,l,e;
1795	list re;
1796	//---------------------------------------------------------------------------
1797	// We need to know which minor corresponds to which variable and to which
1798	// generator of W - therefore we cannot use the function minor()!
1799	//---------------------------------------------------------------------------
1800	//--- choose the generators which we want to differentiate
1801	for(i=1;i<=size(LW);i++)
1802	{
1803	IW=0;
1804	for(l=1;l<=a;l++)
1805	{
1806	if(LW[i][l]!=0){IW[size(IW)+1]=W[l];}
1807	}
1808	//--- choose the variables by which to differentiate and apply diff
1809	for(j=1;j<=size(LV);j++)
1810	{
1811	IV=0;v=0;w=0;
1812	for(l=1;l<=n;l++)
1813	{
1814	if(LV[j][l]!=0)
1815	{
1816	v[size(v)+1]=l;
1817	IV[size(IV)+1]=var(l);
1818	}
1819	else
1820	{
1821	w[size(w)+1]=l;
1822
1823	}
1824	}
1825	A=diff(IV,IW); // appropriate submatrix of Jacobian
1826	//--- if the minor is non-zero, then it might be the one we need
1827	//--- ==> put it in the list of candidates
1828	if(det(A)!=0)
1829	{
1830	v=v[2..size(v)]; // first entry is zero for technical reasons
1831	w=w[2..size(w)]; // first entry is zero for technical reasons
1832	L=inverse_L(A);
1833	L[3]=IW;
1834	L[4]=v;
1835	L[5]=w;
1836	re[size(re)+1]=L;
1837	}
1838	}
1839	}
1840	//---------------------------------------------------------------------------
1841	// return the result
1842	//---------------------------------------------------------------------------
1843	return(re);
1844	}
1845
1846	/////////////////////////////////////////////////////////////////////////////
1847	static proc indexSet(int a, int b)
1848	"Internal procedure - no help and no example available
1849	"
1850	{
1851	//---------------------------------------------------------------------------
1852	// Find all tuples of size b containing pairwise distict elements from a
1853	// list of a elements
1854	//---------------------------------------------------------------------------
1855	//**************************************************************************/
1856	// Combinatorics is expensive in an interpreted language
1857	// ==> shift it into the kernel
1858	//**************************************************************************/
1859	def R=basering;
1860	list L;
1861	ring S=2,x(1..a),dp;
1862	ideal I=maxideal(b);
1863	int i;
1864	ideal J=x(1)^2;
1865	for(i=2;i<=a;i++){J=J,x(i)^2;}
1866	attrib(J,"isSB",1);
1867	I=reduce(I,J);
1868	I=simplify(I,2);
1869	for(i=1;i<=size(I);i++){L[i]=leadexp(I[i]);}
1870	setring R;
1871	return(L);
1872	}
1873
1874	/////////////////////////////////////////////////////////////////////////////
1875
1876	proc DeltaList(list BO)
1877	"USAGE: DeltaList (BO);
1878	@* BO = basic object, a list: ideal W,
1879	@* ideal J,
1880	@* intvec b,
1881	@* list Ex,
1882	@* ideal ab,
1883	@* intvec v,
1884	@* intvec w,
1885	@* matrix M
1886	ASSUME: R = basering, a polynomial ring, W an ideal of R,
1887	@* J = ideal containing W
1888	COMPUTE: Delta-operator iteratively applied to J in the notation of
1889	[Bravo,Encinas,Villamayor]
1890	RETURN: list l of length ((max w-ord) * b),
1891	l[i+1]=Delta^i(J)
1892	EXAMPLE: example DeltaList; shows an example
1893	"
1894	{
1895	//----------------------------------------------------------------------------
1896	// Iteratively apply proc Delta
1897	//----------------------------------------------------------------------------
1898	int i;
1899	list L;
1900	ideal C=BO[2];
1901	while(deg(C[1])!=0)
1902	{
1903	L[size(L)+1]=C;
1904	C=Delta(BO);
1905	BO[2]=C;
1906	}
1907	return(L);
1908	}
1909	example
1910	{
1911	"EXAMPLE:";
1912	echo = 2;
1913	ring R=0,(x,y,z),dp;
1914
1915	ideal W=z^2-x;
1916	ideal J=x*y^2+x^3;
1917	intvec b=1;
1918	list E;
1919	ideal abb=maxideal(1);
1920	intvec v;
1921	intvec w=-1;
1922	matrix M;
1923
1924	list BO=W,J,b,E,abb,v,w,M;
1925
1926	DeltaList(BO);
1927	}
1928	/////////////////////////////////////////////////////////////////////////////
1929	proc CenterBM(list BM)
1930	"USAGE: CenterBM(BM);
1931	@* BM = object related to a presentation,
1932	a list: ideal W,
1933	@* ideal J,
1934	@* intvec b,
1935	@* list Ex,
1936	@* ideal ab,
1937	@* intvec v,
1938	@* intvec w,
1939	@* matrix M
1940	ASSUME: R = basering, a polynomial ring, W an ideal of R,
1941	@* J = ideal containing W
1942	COMPUTE: the center of the next blow-up of BM in the resolution algorithm
1943	of [Bierstone, Milman]
1944	RETURN: list l,
1945	l[1]: ideal describing the center
1946	l[2]: intvec w obtained in the process of determining l[1]
1947	l[3]: intvec b obtained in the process of determining l[1]
1948	l[4]: intmat invmat obtained in the process of determining l[1]
1949	EXAMPLE: example CenterBM; shows an example
1950	"
1951	{
1952	//!!! NOCH NICHT IN BETRIEB
1953	ERROR("Not implemented yet");
1954	int i,j;
1955	intmat tmat[2][1]=0,-1;
1956	//--- re=center,E^- indices, b vector, n vector
1957	list re=ideal(1),BM[7],BM[3],tmat;
1958	ideal J=BM[2];
1959	if(size(J)==0)
1960	{
1961	re[1]=ideal(0);
1962	return(re);
1963	}
1964	//--- find Delta^(b-1)(J)
1965	if(size(reduce(J,std(BM[1])))!=0)
1966	{
1967	list L=DeltaList(BM);
1968	}
1969	else
1970	{
1971	list L;
1972	L[1]=J;
1973	}
1974	if(!defined(debugCenter))
1975	{
1976	int debugCenter;
1977	}
1978	if(debugCenter)
1979	{
1980	"----> In Center: after DeltaList";
1981	"W";
1982	BM[1];
1983	"J";
1984	BM[2];
1985	"The Delta List:";
1986	L;
1987	}
1988	int b=size(L);
1989	if(b==0)
1990	{
1991	//--- if J=W, we do not need to do anything
1992	//--- returning center=1 marks this chart as completed
1993	return(re);
1994	}
1995	//---------------------------------------------------------------------------
1996	// check whether max ord is constant
1997	//---------------------------------------------------------------------------
1998	if((BM[9][2,1]<0)\|\|(BM[9][1,1]>b))
1999	{
2000	//--- we are either at the beginning or the invariant has dropped
2001	intvec tempvec=size(BM[4]);
2002	BM[7]=tempvec;
2003	//!!!! nur fuer hyperflaechen!!!!!!!!
2004	tempvec=b;
2005	BM[3]=tempvec;
2006	//!!!! Ende !!!!!!!!!!!!
2007	kill tempvec;
2008	BM[9][1,1]=b;
2009	BM[9][2,1]=1;
2010	}
2011	//---------------------------------------------------------------------------
2012	// prepare for intersection with E_i
2013	//---------------------------------------------------------------------------
2014	ideal C=L[b];
2015	re[2]=BM[7];
2016	re[3]=BM[3];
2017	BM[2]=C;
2018	if(debugCenter)
2019	{
2020	"----> In Center: before intersection with E_i:";
2021	"bmax:",b;
2022	"Sing(J,bmax):";
2023	C;
2024	"E:";
2025	BO[4];
2026	"list marking a priori known intersection properties:",BO[6];
2027	"index of last element of E^- in E:",BO[7][1];
2028	}
2029	list E=inters_E(BM);
2030	// !!!!!!!!! Drop Redundant fehlt noch!!!!!!!!
2031	//---------------------------------------------------------------------------
2032	// Check whether it is a single point
2033	//---------------------------------------------------------------------------
2034	ideal C1=std(ideal(L[b])+E[1]);
2035	if(dim(C1)==0)
2036	{
2037	if(size(E[4])==1)
2038	{
2039	tmat[1,1]=BM[9][1,1];
2040	tmat[2,1]=BM[9][2,1];
2041	re[4]=tmat;
2042	re[1]=radical(C1);
2043	return(re);
2044	}
2045	}
2046	if(size(BM[9])>2)
2047	{
2048	BM[9][1,2]=E[2];
2049	BM[9][2,2]=1;
2050	}
2051	else
2052	{
2053	intmat tempInt[2][1]=E[2],1;
2054	BM[9]=concatInt(BM[9],tempInt);
2055	kill tempInt;
2056	}
2057	BM[2]=J;
2058	list BM1=dropDim(BM);
2059	list BMlist,hilfList;
2060	ideal hilf;
2061	intvec tempvec;
2062
2063	if(!attrib(BM1[1],"isSB")){BM1[1]=std(BM1[1]);}
2064	for(i=1;i<=size(BM1[4]);i++)
2065	{
2066	hilf=simplify(reduce(BM1[4][i],BM1[1]),2);
2067	if(size(hilf)>1){"Problem with BM1[4]in CenterBM";~;}
2068	hilfList[i]=hilf[1];
2069	}
2070	for(i=1;i<=size(E[4]);i++)
2071	{
2072	BMlist[i]=BM1;
2073	for(j=size(E[5][i]);j>=1;j--)
2074	{
2075	if(E[5][i][j]!=0)
2076	{
2077	BMlist[i][2][size(BMlist[i][2])+1]=hilfList[j];
2078	BMlist[i][3][size(BMlist[i][3])+1]=1;
2079	}
2080	BMlist[i][4]=delete(BMlist[i][4],j);
2081	BMlist[i][6]=deleteInt(BMlist[i][6],j,0);
2082	}
2083	BMlist[i][7]=deleteInt(BMlist[i][7],1,-1);
2084	if(size(BMlist[i][9])>4)
2085	{
2086	intmat tempInt[2][ncols(BMlist[i][9])]=BMlist[i][9];
2087	intmat tempInt2[2][ncols(BMlist[i][9])-2]=
2088	tempInt[1..2,3..ncols(BMlist[i][9])];
2089	BMlist[i][9]=tempInt2;
2090	kill tempInt,tempInt2;
2091	}
2092	else
2093	{
2094	BMlist[i][9]=tmat;
2095	}
2096	}
2097	kill hilfList;list hilfList;
2098	hilfList[1]=CenterTail(BMlist[1],C);
2099	intmat maxmat=hilfList[1][9];
2100	intvec maxiv=E[4][1];
2101	int pos=1;
2102	for(i=2;i<=size(E[4]);i++)
2103	{
2104	hilfList[i]=CenterTail(BMlist[i],C);
2105	if(invGreater(hilfList[i][4]),maxmat,E[4][i],maxiv)
2106	{
2107	maxmat=hilfList[i][4];
2108	maxiv=E[4][i];
2109	pos=i;
2110	}
2111	}
2112	re[1]=hilfList[pos][1];
2113	intmat tempint=BM[9];
2114	intmat tempint1[2][2]=tempint[1..2,1..2];
2115	re[4]=concatInt(tempint1,maxmat);
2116	re[2]=re[2][1],hilfList[pos][2];
2117	re[3]=b;
2118	~;
2119	return(re);
2120	}
2121	example
2122	{ "EXAMPLE:";
2123	echo = 2;
2124	ring R=0,(x,y),dp;
2125
2126	ideal W;
2127	ideal J=x2-y3;
2128	intvec b=1;
2129	list E;
2130	ideal abb=maxideal(1);
2131	intvec v;
2132	intvec w=-1;
2133	matrix M;
2134	intmat invmat[2][1]=0,-1;
2135
2136	list BM=W,J,b,E,abb,v,w,M,invmat;
2137
2138	CenterBM(BM);
2139	}
2140
2141	/////////////////////////////////////////////////////////////////////////////
2142	static proc invGreater(intmat M1, intmat M2, intvec iv1, intvec iv2)
2143	{
2144	// Auxilliary procedure, BM-algorithm
2145	int i;
2146	for(i=1;i<=min(ncols(M1),ncols(M2));i++)
2147	{
2148	if(M1[2,i]==-1)
2149	{
2150	if(M1[1,i]==0){ERROR("Invariant not set");}
2151	if(M2[2,i]!=-1){return(1);}
2152	if(M2[1,i]==0){ERROR("Invariant not set");}
2153	break;
2154	}
2155	else
2156	{
2157	if(M2[2,i]==-1)
2158	{
2159	if(M2[1,i]==0){ERROR("Invariant not set");}
2160	return(0);
2161	}
2162	if(M1[1,i]M2[2,i]!= M2[1,i]M1[2,i])
2163	{
2164	return(M1[1,i]M2[2,i]> M2[1,i]M1[2,i]);
2165	}
2166	}
2167	}
2168	return(iv1>iv2);
2169	}
2170	/////////////////////////////////////////////////////////////////////////////
2171	proc CenterTail(list BM, ideal C)
2172	{
2173	//!!! Auxilliary procedure, BM-algorithm
2174	//!!!!!!!!Rueckgabe im Zentrumsformat
2175	int i,j,bmin;
2176	int alpha=lcm(BM[3]);
2177	vector w;
2178	list re;
2179	if(size(BM[2])==0)
2180	{
2181	re[1]=C+BM[1];
2182	intvec tvec;
2183	re[3]=tvec;
2184	intmat tmat[2][1]=-1,-1;
2185	re[4]=tmat;
2186	tvec=size(BM[4]);
2187	re[2]=tvec;
2188	return(re);
2189	}
2190	for(i=1;i<=size(BM[3]);i++)
2191	{
2192	if(BM[2][i]!=0)
2193	{
2194	w[size(w)+1]=BM[2][i]^(alpha/BM[3][i]);
2195	}
2196	}
2197	module M=w;
2198	intvec satex;
2199	list satList;
2200	for(i=1;i<=size(BM[4]);i++)
2201	{
2202	satList=sat(M,BM[4][i]);
2203	satex[i]=satList[2];
2204	M=satList[1];
2205	}
2206	//!!!!Hilfsobjekt G bilden!!!!!!!!!!!!!!!!
2207	//!!!!Hilfsobjekt H,codim -1 bilden!!!!!!!!!!!!!!!!
2208	//!!!! ???? an welcher stelle?????????
2209	list deltaL;
2210	list BMtemp=BM;
2211	for(i=1;i<=nrows(M[1]);i++)
2212	{
2213	BMtemp[2]=ideal(M[1][i])+BMtemp[1];
2214	deltaL[i]=DeltaList(BMtemp);
2215	if(i==1)
2216	{
2217	bmin=size(deltaL[i]);
2218	}
2219	else
2220	{
2221	bmin=min(size(deltaL[i]),bmin);
2222	}
2223	}
2224	if(bmin==0)
2225	{
2226	re[1]=C+BM[1];
2227	intvec tvec;
2228	re[3]=tvec;
2229	if((BM[9][2,1]==-1)\|\|(BM[9][1,1]!=0))
2230	{
2231	tvec=size(BM[4]);
2232	re[2]=tvec;
2233	intmat tmat[2][1]=0,1;
2234	re[4]=tmat;
2235	}
2236	else
2237	{
2238	re[2]=BM[7];
2239	re[4]=BM[9];
2240	}
2241	return(re);
2242	}
2243	ideal Ctemp=ideal(1);
2244	while(deg(Ctemp[1])==0)
2245	{
2246	Ctemp=C;
2247	for(i=1;i<=nrows(M[1]);i++)
2248	{
2249	Ctemp=Ctemp,deltaL[i][bmin];
2250	}
2251	Ctemp=std(Ctemp);
2252	bmin--;
2253	if(bmin==0){ERROR("empty set");}
2254	}
2255	//!!!!!!!!!!!!!Invariante ist bmin/alpha
2256	// naechster Eintrag s_i wie in CenterBM
2257	// dann dropDim ....
2258	}
2259
2260	/////////////////////////////////////////////////////////////////////////////
2261	static proc deleteInt(intvec v,int i,int ini)
2262	{
2263	//!!! Should be in kernel of Singular
2264	//--- delete i-th entry in intvec v,
2265	//--- if necessary reinitializing v with value ini
2266	int s=size(v);
2267	intvec w;
2268	if((i<s)&&(i>1)){w=v[1..i-1],v[i+1..s];}
2269	if(s==1){w=ini;return(w);}
2270	if(i==1){w=v[2..s];}
2271	if(i==s){w=v[1..s-1];}
2272	return(w);
2273	}
2274	/////////////////////////////////////////////////////////////////////////////
2275	static proc concatInt(intmat A, intmat B)
2276	{
2277	//!!! Should be in kernel of Singular
2278	//--- concatenate two intmats
2279	if(nrows(A)!=nrows(B)){ERROR("could not concat, wrong number of rows");}
2280	intmat tempmat[nrows(A)][ncols(A)+ncols(B)];
2281	tempmat[1..nrows(A),1..ncols(A)]=A[1..nrows(A),1..ncols(A)];
2282	tempmat[1..nrows(A),ncols(A)+1..ncols(tempmat)]=B[1..nrows(A),1..ncols(B)];
2283	return(tempmat);
2284	}
2285	/////////////////////////////////////////////////////////////////////////////
2286	proc dropDim(list BM)
2287	{
2288	ERROR("Not implemented yet");
2289	}
2290	/////////////////////////////////////////////////////////////////////////////
2291	proc CenterBO(list BO,list #)
2292	"USAGE: CenterBO(BO);
2293	@* BO = basic object, a list: ideal W,
2294	@* ideal J,
2295	@* intvec b,
2296	@* list Ex,
2297	@* ideal ab,
2298	@* intvec v,
2299	@* intvec w,
2300	@* matrix M
2301	ASSUME: R = basering, a polynomial ring, W an ideal of R,
2302	@* J = ideal containing W
2303	COMPUTE: the center of the next blow-up of BO in the resolution algorithm
2304	of [Bravo,Encinas,Villamayor]
2305	RETURN: list l,
2306	l[1]: ideal describing the center@*
2307	l[2]: intvec w obtained in the process of determining l[1]@*
2308	l[3]: intvec b obtained in the process of determining l[1]@*
2309	l[4]: intvec inv obtained in the process of determining l[1]
2310	EXAMPLE: example CenterBO; shows an example
2311	"
2312	{
2313	//---------------------------------------------------------------------------
2314	// Initialization and sanity checks
2315	//---------------------------------------------------------------------------
2316	int i,bo7save;
2317	intvec tvec;
2318	//--- re=center,E^- indices, b vector, n vector
2319	list re=ideal(1),BO[7],BO[3],tvec;
2320	ideal J=BO[2];
2321	if(size(J)==0)
2322	{
2323	re[1]=ideal(0);
2324	return(re);
2325	}
2326	//--- find Delta^(b-1)(J)
2327	if(size(reduce(J,std(BO[1])))!=0)
2328	{
2329	list L=DeltaList(BO);
2330	}
2331	else
2332	{
2333	list L;
2334	L[1]=J;
2335	}
2336	if(!defined(debugCenter))
2337	{
2338	int debugCenter;
2339	}
2340	if(debugCenter)
2341	{
2342	"----> In Center: after DeltaList";
2343	"W";
2344	BO[1];
2345	"J";
2346	BO[2];
2347	"The Delta List:";
2348	L;
2349	}
2350	int b=size(L);
2351	if(b==0)
2352	{
2353	//--- if J=W, we do not need to do anything
2354	//--- returning center=1 marks this chart as completed
2355	return(re);
2356	}
2357	//---------------------------------------------------------------------------
2358	// check whether max w-ord is constant
2359	//---------------------------------------------------------------------------
2360	if(b==BO[3][1])
2361	{
2362	//--- max w-ord is constant
2363	if(BO[7][1]==-1)
2364	{
2365	//--- first got its value in the previous step ==> initialize BO[7]
2366	tvec[1]=size(BO[4])-1;
2367	for(i=2;i<=size(BO[7]);i++)
2368	{
2369	tvec[i]=BO[7][i];
2370	}
2371	re[2]=tvec;
2372	BO[7]=tvec;
2373	}
2374	}
2375	else
2376	{
2377	//--- max w-ord changed ==> reset BO[7], correct BO[3]
2378	tvec[1]=-1;
2379	re[2]=tvec;
2380	BO[7]=tvec;
2381	tvec[1]=b;
2382	BO[3]=tvec;
2383	if(defined(invSat))
2384	{
2385	invSat[2]=intvec(0);
2386	}
2387	}
2388	re[3]=BO[3];
2389	//---------------------------------------------------------------------------
2390	// reduce from case 2 to case 1 of [Bravo, Encinas, Villamayor]
2391	//---------------------------------------------------------------------------
2392	ideal C=L[b];
2393	BO[2]=C;
2394	if(debugCenter)
2395	{
2396	"----> In Center: before intersection with E_i:";
2397	"bmax:",b;
2398	"Sing(J,bmax):";
2399	C;
2400	"E:";
2401	BO[4];
2402	"list marking a priori known intersection properties:",BO[6];
2403	"index of last element of E^- in E:",BO[7][1];
2404	}
2405	//--- is intermediate result in iteration already good?
2406	//--- return it to calling proc CenterBO
2407	if(size(#)>0)
2408	{
2409	if(#[1]==2)
2410	{
2411	re[1]=C;
2412	kill tvec;
2413	intvec tvec=re[2][1];
2414	re[2]=tvec;
2415	kill tvec;
2416	intvec tvec=re[3][1];
2417	re[3]=tvec;
2418	kill tvec;
2419	intvec tvec;
2420	re[4]=tvec;
2421	return(re);
2422	}
2423	}
2424	//--- do the reduction to case 1
2425	list E=inters_E(BO);
2426	//--- if J is smooth and not too many E_i intersect simultaneously, let us
2427	//--- try to drop redundant components of the candidate for the center
2428	if((b==1)&&(size(E)>3))
2429	{
2430	//--- if J is not smooth we do not want to drop any information
2431	if((size(E[4])>0) && (dim(std(slocusE(BO[2])))<0))
2432	{
2433	//--- BO[2]==J because b==1
2434	//--- DropRedundant is the counterpart to DropCoeff
2435	//--- do not leave out one of them separately!!!
2436	E=DropRedundant(BO,E);
2437	if(size(E)==1)
2438	{
2439	kill tvec;
2440	intvec tvec=re[2][1];
2441	re[2]=tvec;
2442	tvec[1]=re[3][1];
2443	re[3]=tvec;
2444	tvec[1]=re[4][1];
2445	re[4]=tvec;
2446	re[1]=E[1];
2447	return(re);
2448	}
2449	}
2450	}
2451	//--- set n correctly
2452	if(E[2]<BO[9][1])
2453	{
2454	//--- if n dropped, the subsequent Coeff-object will not be the same
2455	//--- ===> set BO[3][2] to zero to make sure that no previous data is used
2456	if(defined(tvec)) {kill tvec;}
2457	intvec tvec=BO[7][1],-1;
2458	BO[7]=tvec;
2459	tvec=BO[3][1],0;
2460	BO[3]=tvec;
2461	if(defined(invSat))
2462	{
2463	invSat[2]=intvec(0);
2464	}
2465	}
2466	re[4][1]=E[2];
2467	C=E[1]^b+J;
2468	C=mstd(C)[2];
2469	ideal C1=std(ideal(L[b])+E[1]);
2470	if(debugCenter)
2471	{
2472	"----> In Center: reduction of case 2 to case 1";
2473	"Output of inters_E, after dropping redundant components:";
2474	E;
2475	"result of intersection with E^-, i.e.(E^-)^b+J:";
2476	C;
2477	"candidate for center:";
2478	C1;
2479	}
2480	//---------------------------------------------------------------------------
2481	// Check whether we have a hypersurface component
2482	//---------------------------------------------------------------------------
2483	if(dim(C1)==dim(std(BO[1]))-1)
2484	{
2485	if((size(reduce(J,C1))==0)&&(size(reduce(C1,std(J)))==0))
2486	{
2487	//--- C1 equals J and is of codimension 1 in W
2488	re[1]=C1;
2489	}
2490	else
2491	{
2492	//--- C1 has a codimension 1 (in W) component
2493	re[1]=std(equiRadical(C1));
2494	}
2495	kill tvec;
2496	intvec tvec=re[2][1];
2497	re[2]=tvec;
2498	tvec[1]=re[3][1];
2499	re[3]=tvec;
2500	tvec[1]=re[4][1];
2501	re[4]=tvec;
2502	//--- is the codimension 1 component a good choice or do we need to reset
2503	//--- the information from the previous steps
2504	if(transversalT(re[1],BO[4]))
2505	{
2506	if(size(E)>2)
2507	{
2508	if(E[3]>E[2])
2509	{
2510	if(defined(shortcut)){kill shortcut;}
2511	list shortcut=ideal(0),size(BO[4]),BO[7];
2512	export(shortcut);
2513	}
2514	}
2515	return(re);
2516	}
2517
2518	ERROR("reset in Center, please send the example to the authors.");
2519	}
2520	//---------------------------------------------------------------------------
2521	// Check whether it is a single point
2522	//---------------------------------------------------------------------------
2523	if(dim(C1)==0)
2524	{
2525	C1=std(radical(C1));
2526	if(vdim(C1)==1)
2527	{
2528	//--- C1 is one point
2529	re[1]=C1;
2530	kill tvec;
2531	intvec tvec=re[2][1];
2532	re[2]=tvec;
2533	kill tvec;
2534	intvec tvec=re[3][1];
2535	re[3]=tvec;
2536	return(re);
2537	}
2538	}
2539	//---------------------------------------------------------------------------
2540	// Prepare input for forming the Coeff-Ideal
2541	//---------------------------------------------------------------------------
2542	BO[2]=C;
2543	if(size(BO[2])>5)
2544	{
2545	BO[2]=mstd(BO[2])[2];
2546	}
2547	//--- drop leading entry of BO[3]
2548	tvec=BO[3];
2549	if(size(tvec)>1)
2550	{
2551	tvec=tvec[2..size(tvec)];
2552	BO[3]=tvec;
2553	}
2554	else
2555	{
2556	BO[3][1]=0;
2557	}
2558	tvec=BO[9];
2559	if(size(tvec)>1)
2560	{
2561	tvec=tvec[2..size(tvec)];
2562	BO[9]=tvec;
2563	}
2564	else
2565	{
2566	BO[9][1]=0;
2567	}
2568	bo7save=BO[7][1]; // original value needed for result
2569	if(defined(shortcut))
2570	{
2571	if((bo7save!=shortcut[3][1])&&(size(shortcut[3])!=1))
2572	{
2573	kill shortcut;
2574	}
2575	else
2576	{
2577	shortcut[2]=shortcut[2]-bo7save;
2578	tvec=shortcut[3];
2579	if(size(tvec)>1)
2580	{
2581	tvec=tvec[2..size(tvec)];
2582	shortcut[3]=tvec;
2583	}
2584	else
2585	{
2586	shortcut[3]=intvec(shortcut[2]);
2587	}
2588	}
2589	}
2590	if(BO[7][1]>-1)
2591	{
2592	//--- drop E^- and the corresponding information from BO[6]
2593	for(i=1;i<=BO[7][1];i++)
2594	{
2595	BO[4]=delete(BO[4],1);
2596	intvec bla1=BO[6];
2597	BO[6]=intvec(bla1[2..size(bla1)]);
2598	kill bla1;
2599	}
2600	//--- drop leading entry of BO[7]
2601	tvec=BO[7];
2602	if(size(tvec)>1)
2603	{
2604	tvec=tvec[2..size(tvec)];
2605	BO[7]=tvec;
2606	}
2607	else
2608	{
2609	BO[7][1]=-1;
2610	}
2611	}
2612	else
2613	{
2614	if(BO[7][1]==-1)
2615	{
2616	list tplist;
2617	BO[4]=tplist;
2618	kill tplist;
2619	}
2620	}
2621	if(debugCenter)
2622	{
2623	"----> In Center: Input to Coeff";
2624	"b:",b;
2625	"BO:";
2626	BO;
2627	}
2628	//--- prepare the third entry of the invariant tuple
2629	int invSatSave=invSat[2][1];
2630	tvec=invSat[2];
2631	if(size(tvec)>1)
2632	{
2633	tvec=tvec[2..size(tvec)];
2634	invSat[2]=tvec;
2635	}
2636	else
2637	{
2638	invSat[2][1]=0;
2639	}
2640	//---------------------------------------------------------------------------
2641	// Form the Coeff-ideal, if possible and useful; otherwise use the previous
2642	// candidate for the center
2643	//---------------------------------------------------------------------------
2644	list BO1=Coeff(BO,b);
2645	if(debugCenter)
2646	{
2647	"----> In Center: Output of Coeff";
2648	BO1;
2649	}
2650	//--- Coeff returns int if something went wrong
2651	if(typeof(BO1[1])=="int")
2652	{
2653	if(BO1[1]==0)
2654	{
2655	//--- Coeff ideal was already resolved
2656	re[1]=C1;
2657	return(re);
2658	}
2659	else
2660	{
2661	//--- no global hypersurface found
2662	re=CoverCenter(BO,b,BO1[2]);
2663	kill tvec;
2664	intvec tvec=invSatSave;
2665	for(i=1;i<=size(invSat[2]);i++)
2666	{
2667	tvec[i+1]=invSat[2][i];
2668	}
2669	invSat[2]=tvec;
2670	return(re);
2671	}
2672	}
2673	int coeff_invar;
2674	ideal Idropped=1;
2675	//--- if b=1 drop redundant components of the Coeff-ideal
2676	if(b==1)
2677	{
2678	//--- Counterpart to DropRedundant -- do not leave out one of them separately
2679	Idropped=DropCoeff(BO1); // blow-up in these components
2680	// is unnecessary
2681	}
2682	//--- to switch off DropCoeff, set Idropped=1;
2683	BO1[2]=sat(BO1[2],Idropped)[1];
2684	if(deg(BO1[2][1])==0)
2685	{
2686	//--- Coeff ideal is trivial
2687	C1=radical(C1);
2688	ideal C2=sat(C1,Idropped)[1];
2689	if(deg(std(C2)[1])!=0)
2690	{
2691	C1=C2;
2692	}
2693	//Aenderung: Strategie: nur im Notfall ganze except. Divisoren
2694	if(deg(std(BO1[2])[1])==0)
2695	{
2696	list BOtemp=BO;
2697	int bo17save=BO1[7][1];
2698	BOtemp[7]=0;
2699	BO1=Coeff(BOtemp,b,int(0));
2700	BO1[2]=sat(BO1[2],Idropped)[1];
2701	if(deg(std(BO1[2])[1])==0)
2702	{
2703	//--- there is really nothing left to do for the Coeff ideal
2704	//--- the whole original BO1[2], i.e. Idropped, is the upcoming center
2705	re[1]=Idropped;
2706	re[2]=intvec(bo7save);
2707	re[3]=intvec(b);
2708	re[4]=intvec(E[2]);
2709	return(re);
2710	}
2711	if(deg(std(slocus(radical(BO1[2])))[1])==0)
2712	{
2713	re[1]=BO1[2];
2714	// re[2]=intvec(bo7save,BO1[7][1]);
2715	re[2]=intvec(bo7save,bo17save);
2716	re[3]=intvec(b,1);
2717	re[4]=intvec(E[2],1);
2718	invSat[2]=intvec(1,0);
2719	return(re);
2720	}
2721	//!!! effizienter machen???
2722	list pr=primdecGTZ(BO1[2]);
2723	ideal Itemp1=1;
2724	int aa,bb;
2725	for(aa=1;aa<=size(pr);aa++)
2726	{
2727	if(dim(std(pr[aa][2])) < (dim(std(BO1[1]))-1))
2728	{
2729	//--- drop components which are themselves exceptional diviosrs
2730	Itemp1=intersect(Itemp1,pr[aa][1]);
2731	}
2732	}
2733	if(deg(std(Itemp1)[1])!=0)
2734	{
2735	//--- treat the remaining components of the weak Coeff ideal
2736	BO1[2]=Itemp1;
2737	}
2738	BO1[7]=BO[7];
2739	for(aa=1;aa<=size(BO1[4]);aa++)
2740	{
2741	if(deg(std(BO1[4][aa])[1])==0){aa++;continue;}
2742	if(defined(satlist)){kill satlist;}
2743	list satlist=sat(BO1[2],BO1[4][aa]+BO1[1]);
2744	if(deg(std(satlist[1])[1])==0)
2745	{
2746	coeff_invar++;
2747	if(satlist[2]!=0)
2748	{
2749	for(bb=1;bb<=satlist[2]-1;bb++)
2750	{
2751	BO1[2]=quotient(BO1[2],BO1[4][aa]+BO1[1]);
2752	}
2753	}
2754	else
2755	{
2756	ERROR("J of temporary object had unexpected value;
2757	please send this example to the authors.");
2758	}
2759	}
2760	else
2761	{
2762	BO1[2]=satlist[1];
2763	}
2764	}
2765	if(deg(std(Itemp1)[1])==0)
2766	{
2767	re[1]=BO1[2];
2768	re[2]=intvec(bo7save,BO1[7][1]);
2769	re[3]=intvec(b,1);
2770	re[4]=intvec(E[2],1);
2771	invSat[2]=intvec(1,0);
2772	return(re);
2773	}
2774	kill aa,bb;
2775	}
2776	}
2777	if(invSatSave<coeff_invar)
2778	{
2779	invSatSave=coeff_invar;
2780	}
2781	//---------------------------------------------------------------------------
2782	// Determine Center of Coeff-ideal and use it as the new center
2783	//---------------------------------------------------------------------------
2784	if(!defined(templist))
2785	{
2786	if(size(BO1[2])>5)
2787	{
2788	BO1[2]=mstd(BO1[2])[2];
2789	}
2790	list templist=CenterBO(BO1,2);
2791	//--- only a sophisticated guess of a good center computed by
2792	//--- leaving center before intersection with the E_i.
2793	//--- whether the guess was good, is stored in 'good'.
2794	//--- (this variant saves charts in cases like the Whitney umbrella)
2795	list E0,E1;
2796
2797	ideal Cstd=std(radical(templist[1]));
2798	int good=((deg(std(slocusE(Cstd))[1])==0)&&(dim(std(BO1[2]))<=2));
2799	//if(defined(satlist)){good=0;}
2800	if(good)
2801	{
2802	for(i=1;i<=size(BO[4]);i++)
2803	{
2804	if((deg(BO[4][i][1])>0)&&(size(reduce(BO[4][i],Cstd))!=0))
2805	{
2806	E0[size(E0)+1]=BO[4][i]+Cstd;
2807	E1[size(E1)+1]=BO[4][i];
2808	}
2809	}
2810	good=transversalT(Cstd,E1);
2811	if(good)
2812	{
2813	good=normalCross(E0);
2814	}
2815	}
2816	if(good)
2817	{
2818	list templist2=CenterBO(BO1,1);
2819	if(dim(std(templist2[1]))!=dim(Cstd))
2820	{
2821	templist[1]=Cstd;
2822	if(defined(shortcut)){kill shortcut;}
2823	list shortcut=ideal(0),size(BO1[4]),templist[2];
2824	export(shortcut);
2825	}
2826	else
2827	{
2828	templist=templist2;
2829	}
2830	kill templist2;
2831	}
2832	else
2833	{
2834	//--- sophisticated guess was wrong, follow Villamayor's approach
2835	kill templist;
2836	list templist=CenterBO(BO1,1);
2837	}
2838	}
2839	if((dim(std(templist[1]))==dim(std(BO1[1]))-1)
2840	&&(size(templist[4])==1))
2841	{
2842	if(templist[4][1]==0)
2843	{
2844	for(i=1;i<=size(BO1[4]);i++)
2845	{
2846	if(size(reduce(templist[1],std(BO1[4][i])))==0)
2847	{
2848	templist[4][1]=1;
2849	break;
2850	}
2851	}
2852	}
2853	}
2854	//!!! subsequent line should be deleted
2855	//if(defined(satlist)){templist[3][1]=BO[3][1];}
2856	if(debugCenter)
2857	{
2858	"----> In Center: Iterated Center returned:";
2859	templist;
2860	}
2861	//--------------------------------------------------------------------------
2862	// set up the result and return it
2863	//--------------------------------------------------------------------------
2864	re[1]=templist[1];
2865	kill tvec;
2866	intvec tvec;
2867	tvec[1]=bo7save;
2868	for(i=1;i<=size(templist[2]);i++)
2869	{
2870	tvec[i+1]=templist[2][i];
2871	}
2872	re[2]=tvec;
2873	if(defined(shortcut))
2874	{
2875	shortcut[2]=shortcut[2]+bo7save;
2876	shortcut[3]=tvec;
2877	}
2878	kill tvec;
2879	intvec tvec;
2880	tvec[1]=invSatSave;
2881	for(i=1;i<=size(invSat[2]);i++)
2882	{
2883	tvec[i+1]=invSat[2][i];
2884	}
2885	invSat[2]=tvec;
2886	kill tvec;
2887	intvec tvec;
2888	tvec[1]=b;
2889	for(i=1;i<=size(templist[3]);i++)
2890	{
2891	tvec[i+1]=templist[3][i];
2892	}
2893	re[3]=tvec;
2894	kill tvec;
2895	intvec tvec;
2896	tvec[1]=E[2];
2897	for(i=1;i<=size(templist[4]);i++)
2898	{
2899	tvec[i+1]=templist[4][i];
2900	}
2901	re[4]=tvec;
2902
2903	return(re);
2904	}
2905	example
2906	{ "EXAMPLE:";
2907	echo = 2;
2908	ring R=0,(x,y),dp;
2909
2910	ideal W;
2911	ideal J=x2-y3;
2912	intvec b=1;
2913	list E;
2914	ideal abb=maxideal(1);
2915	intvec v;
2916	intvec w=-1;
2917	matrix M;
2918
2919	list BO=W,J,b,E,abb,v,w,M,v;
2920
2921	CenterBO(BO);
2922	}
2923	//////////////////////////////////////////////////////////////////////////////
2924	static proc CoverCenter(list BO,int b, ideal Jb)
2925	{
2926	//----------------------------------------------------------------------------
2927	// Initialization
2928	//----------------------------------------------------------------------------
2929	def R=basering;
2930	int i,j,k;
2931	intvec merk,merk2,maxv,fvec;
2932	list L,ceList,re;
2933	ceList[1]=ideal(0);
2934	poly @p,@f;
2935	ideal K,dummy;
2936	if(!attrib(BO[2],"isSB"))
2937	{
2938	BO[2]=std(BO[2]);
2939	}
2940	for(i=1;i<=size(Jb);i++)
2941	{
2942	list tempmstd=mstd(slocus(Jb[i]));
2943	if(size(tempmstd[1])>size(tempmstd[2]))
2944	{
2945	dummy=tempmstd[2];
2946	}
2947	else
2948	{
2949	dummy=tempmstd[1];
2950	}
2951	kill tempmstd;
2952	L[i]=dummy;
2953	K=K,dummy;
2954	}
2955	K=simplify(K,2);
2956	//---------------------------------------------------------------------------
2957	// The intersection of the singular loci of the L[i] is empty.
2958	// Find a suitable open covering of the affine chart, such that a global
2959	// hypersurface can be found in each open set.
2960	//---------------------------------------------------------------------------
2961	matrix M=lift(K,ideal(1));
2962	j=1;
2963	for(i=1;i<=nrows(M);i++)
2964	{
2965	if(M[i,1]!=0)
2966	{
2967	merk[size(merk)+1]=i;
2968	fvec[size(merk)]=j;
2969	}
2970	if((i-k)==size(L[j]))
2971	{
2972	k=i;
2973	j++;
2974	}
2975	}
2976	//--------------------------------------------------------------------------
2977	// Find a candidate for the center in each open set
2978	//--------------------------------------------------------------------------
2979	//--- first entry of merk is 0 by construction of merk
2980	for(i=2;i<=size(merk);i++)
2981	{
2982	//--- open set is D(@p)
2983	@p=K[merk[i]];
2984	//--- hypersurface is V(@f)
2985	@f=Jb[fvec[i]];
2986	execute("ring R1=("+charstr(R)+"),(@y,"+varstr(R)+"),dp;");
2987	poly p=imap(R,@p);
2988	poly f=imap(R,@f);
2989	list @ce;
2990	list BO=imap(R,BO);
2991	BO[1]=BO[1]+ideal(@y*p-1);
2992	BO[2]=BO[2]+ideal(@y*p-1);
2993	for(j=1;j<=size(BO[4]);j++)
2994	{
2995	BO[4][j]=BO[4][j]+ideal(@y*p-1);
2996	}
2997	//--- like usual Coeff, but hypersurface is already known
2998	list BO1=SpecialCoeff(BO,b,f);
2999	//--- special situation in SpecialCoeff are marked by an error code of
3000	//--- type int
3001	if(typeof(BO1[1])=="int")
3002	{
3003	if(BO1[1]==0)
3004	{
3005	//--- Coeff ideal was already resolved
3006	@ce[1]=BO[2];
3007	@ce[2]=BO[7];
3008	@ce[3]=BO[3];
3009	}
3010	else
3011	{
3012	if(BO[3]!=0)
3013	{
3014	//--- intersections with E do not meet conditions ==> reset
3015	ERROR("reset in Coeff, please send the example to the autors");
3016	}
3017	}
3018	}
3019	else
3020	{
3021	//--- now do the recursion as usual
3022	@ce=CenterBO(BO1);
3023	}
3024	//---------------------------------------------------------------------------
3025	// Go back to the whole affine chart and form a suitable union of the
3026	// candidates
3027	//---------------------------------------------------------------------------
3028	//--- pass from open set to the whole affine chart by taking the closure
3029	@ce[1]=eliminate(@ce[1],@y);
3030	setring R;
3031	ceList[i]=imap(R1,@ce);
3032	//--- set up invariant vector and determine maximum value of it
3033	if(size(ceList[i][3])==size(ceList[i][4]))
3034	{
3035	kill merk2,maxv;
3036	intvec merk2,maxv;
3037	for(j=1;j<=size(ceList[i][3]);j++)
3038	{
3039	merk2[2*j-1]=ceList[i][3][j];
3040	merk2[2*j]=ceList[i][4][j];
3041	ceList[i][5]=merk2;
3042	if(maxv<merk2)
3043	{
3044	maxv=merk2;
3045	}
3046	}
3047	}
3048	else
3049	{
3050	ERROR("This situation should not occur, please send the example
3051	to the authors.");
3052	}
3053	kill R1;
3054	}
3055	kill merk2;
3056	intvec merk2=-2;
3057	//--- form the union of the components of the center with maximum invariant
3058	for(i=1;i<=size(ceList);i++)
3059	{
3060	if(size(reduce(ceList[i][1],BO[2]))==0)
3061	{
3062	//--- already resolved ==> ignore
3063	i++;
3064	continue;
3065	}
3066	if(ceList[i][5]==maxv)
3067	{
3068	if(merk2!=ceList[i][2])
3069	{
3070	//--- E^- not of the same size as before resp. initialization
3071	if(merk2[1]==-2)
3072	{
3073	//--- initialization: save size of E^-
3074	merk2=ceList[i][2];
3075	re[1]=ceList[i][1];
3076	re[2]=ceList[i][2];
3077	re[3]=ceList[i][3];
3078	re[4]=ceList[i][4];
3079	}
3080	else
3081	{
3082	//--- otherwise ignore
3083	i++;
3084	continue;
3085	}
3086	}
3087	else
3088	{
3089	re[1]=intersect(re[1],ceList[i][1]);
3090	}
3091	}
3092	}
3093	//--------------------------------------------------------------------------
3094	// Perform last checks and return the result
3095	//--------------------------------------------------------------------------
3096	if(size(re)!=4)
3097	{
3098	//--- oops: already resolved in all open sets
3099	re[1]=BO[2];
3100	re[2]=-1;
3101	re[3]=0;
3102	re[4]=intvec(0);
3103	}
3104	return(re);
3105	}
3106	//////////////////////////////////////////////////////////////////////////////
3107	static proc SpecialCoeff(list BO,int b,poly f)
3108	{
3109	//----------------------------------------------------------------------------
3110	// Coeff with given hypersurface -- no checks of the hypersurface performed
3111	//----------------------------------------------------------------------------
3112	int i,ch;
3113	int e=int(factorial(b));
3114	ideal C;
3115	list L=DeltaList(BO);
3116	int d=size(L);
3117	//--- set up ideal
3118	for(i=0;i<b;i++)
3119	{
3120	C=C+L[i+1]^(e/(b-i));
3121	}
3122	//--- move to hypersurface V(Z)
3123	ideal Z=f;
3124	C=C,Z;
3125	BO[1]=BO[1]+Z;
3126	BO[2]=C;
3127	for(i=1;i<=size(BO[4]);i++)
3128	{
3129	BO[6][i]=0; // reset intersection indicator
3130	BO[4][i]=BO[4][i]+Z; // intersect the E_i
3131	BO[2]=sat(BO[2],BO[4][i]+BO[1])[1];
3132	// "strict transform" of J w.r.t E, not "total"
3133	}
3134	return(BO);
3135	}
3136
3137	//////////////////////////////////////////////////////////////////////////////
3138	static proc DropCoeff(list BO)
3139	"Internal procedure - no help and no example available
3140	"
3141	{
3142	//--- Initialization
3143	int i,j;
3144	list pr=minAssGTZ(BO[2]);
3145	ideal I=BO[2];
3146	ideal Itemp;
3147	ideal Idropped=1;
3148	//--- Tests
3149	for(i=1;i<=size(pr);i++)
3150	{
3151	if(i>size(pr))
3152	{
3153	//--- the continue statement does not test the loop condition sigh
3154	break;
3155	}
3156	if(deg(std(slocus(pr[i]))[1])!=0)
3157	{
3158	//--- this component is singular ===> we still need it
3159	i++;
3160	continue;
3161	}
3162	Itemp=sat(I,pr[i])[1];
3163	if(deg(std(Itemp+pr[i])[1])!=0)
3164	{
3165	//--- this component is not disjoint from the other ones ===> we still need it
3166	i++;
3167	continue;
3168	}
3169	if(!transversalT(pr[i],BO[4]))
3170	{
3171	//--- this component does not meet one of the remaining E_i transversally
3172	//--- ===> we still need it
3173	i++;
3174	continue;
3175	}
3176	if(!normalCross(BO[4],pr[i]))
3177	{
3178	//--- this component is not normal crossing with the remaining E_i
3179	//--- ===> we still need it
3180	i++;
3181	continue;
3182	}
3183	if(defined(EE)){kill EE;}
3184	list EE;
3185	for(j=1;j<=size(BO[4]);j++)
3186	{
3187	EE[j]=BO[4][j]+pr[i];
3188	}
3189	if(!normalCross(EE))
3190	{
3191	//--- we do not have a normal crossing situation for this component after all
3192	//--- ===> we still need it
3193	i++;
3194	continue;
3195	}
3196	Idropped=intersect(Idropped,pr[i]);
3197	I=Itemp;
3198	}
3199	return(Idropped);
3200	}
3201
3202	//////////////////////////////////////////////////////////////////////////////
3203	static proc DropRedundant(list BO,list E)
3204	"Internal procedure - no help and no example available
3205	"
3206	{
3207	//---------------------------------------------------------------------------
3208	// Initialization and sanity checks
3209	//---------------------------------------------------------------------------
3210	int ii,jj,kkdiff,nonnormal,ok;
3211	ideal testid,dummy;
3212	ideal center;
3213	intvec transverse,dontdrop,zerovec;
3214	transverse[size(BO[4])]=0;
3215	dontdrop[size(E[4])]=0;
3216	zerovec[size(E[4])]=0;
3217	ideal J=BO[2];
3218	int dimJ=dim(std(BO[2]));
3219	list templist;
3220	if(size(E)<5)
3221	{
3222	//--- should not occur
3223	return(E);
3224	}
3225	for(ii=1;ii<=BO[7][1];ii++)
3226	{
3227	if(BO[6][ii]==2)
3228	{
3229	kkdiff++;
3230	}
3231	}
3232	int expDim=dimJ-E[2]+kkdiff;
3233	if(size(E)==6)
3234	{
3235	nonnormal=E[6];
3236	}
3237	//---------------------------------------------------------------------------
3238	// if dimJ were smaller than E[2], we would not have more than 3 entries in
3239	// in the list E
3240	// if dimJ is also at least E[3] and nonnormal is 0, we only need to test that
3241	// * the intersection is of the expected dimension
3242	// * the intersections of the BO[4][i] and J are normal crossing
3243	// * the elements of E^+ have no influence (is done below)
3244	//---------------------------------------------------------------------------
3245	if((E[3]<=dimJ)&&(!nonnormal))
3246	{
3247	ideal bla=E[1]+BO[2]+BO[1];
3248	bla=radical(bla);
3249	bla=mstd(bla)[2];
3250
3251	if(dim(std(slocusE(bla)))<0)
3252	{
3253	if(transversalT(J,BO[4]))
3254	{
3255	ok=1;
3256	if(E[2]==E[3])
3257	{
3258	//--- no further intersection with elements from E^+
3259	for(ii=1;ii<=size(E[4]);ii++)
3260	{
3261	if(dim_slocus(BO[2]+E[4][ii])!=-1)
3262	{
3263	dontdrop[ii]=1;
3264	}
3265	}
3266	if(dontdrop==zerovec)
3267	{
3268	list relist;
3269	relist[1]=std(J);
3270	return(relist);
3271	}
3272	}
3273	}
3274	}
3275	}
3276	//---------------------------------------------------------------------------
3277	// now check whether the E_i actually occurring in the intersections meet
3278	// J transversally (one by one) and mark those elements of E[4] where it is
3279	// not the case
3280	//---------------------------------------------------------------------------
3281	if(!ok)
3282	{
3283	for(ii=1;ii<=size(E[5]);ii++)
3284	{
3285	//--- test the ii-th tuple of E[4] resp. its indices E[5]
3286	for(jj=1;jj<=size(E[5][ii]);jj++)
3287	{
3288	//--- if E[5][ii][jj]==1, E_jj is involved in E[4][ii]
3289	if(E[5][ii][jj]==1)
3290	{
3291	//--- transversality not yet determined
3292	if(transverse[jj]==0)
3293	{
3294	templist[1]=BO[4][jj];
3295	if(transversalT(BO[2],templist))
3296	{
3297	transverse[jj]=1;
3298	}
3299	else
3300	{
3301	transverse[jj]=-1;
3302	dontdrop[ii]=1;
3303	}
3304	}
3305	else
3306	{
3307	//--- already computed transversality
3308	if(transverse[jj]<0)
3309	{
3310	dontdrop[ii]=1;
3311	}
3312	}
3313	}
3314	}
3315	}
3316	}
3317	//---------------------------------------------------------------------------
3318	// if one of the non-marked tuples from E^- in E[4] has an intersection
3319	// of the expected dimension and does not meet any E_i from E^+
3320	// - except the ones which are met trivially - , it should be
3321	// dropped from the list.
3322	// it can also be dropped if an intersection occurs and normal crossing has
3323	// been checked.
3324	//---------------------------------------------------------------------------
3325	for(ii=1;ii<=size(E[4]);ii++)
3326	{
3327	//--- if E[4][ii] does not have transversal intersections, we cannot drop it
3328	if(dontdrop[ii]==1)
3329	{
3330	ii++;
3331	continue;
3332	}
3333	//--- testing ii-th tuple from E[4]
3334	testid=BO[1]+BO[2]+E[4][ii];
3335	if(dim(std(testid))!=expDim)
3336	{
3337	//--- not expected dimension
3338	dontdrop[ii]=1;
3339	ii++;
3340	continue;
3341	}
3342	testid=mstd(testid)[2];
3343
3344	if(dim(std(slocusE(testid)))>=0)
3345	{
3346	//--- not smooth, i.e. more than one component which intersect
3347	dontdrop[ii]=1;
3348	ii++;
3349	continue;
3350	}
3351	//--- if E^+ is empty, we are done; otherwise check intersections with E^+
3352	if(BO[7][1]!=-1)
3353	{
3354	if(defined(pluslist)){kill pluslist;}
3355	list pluslist;
3356	for(jj=BO[7][1]+1;jj<=size(BO[4]);jj++)
3357	{
3358	dummy=BO[4][jj]+testid;
3359	dummy=std(dummy);
3360	if(expDim==dim(dummy))
3361	{
3362	//--- intersection has wrong dimension
3363	dontdrop[ii]=1;
3364	break;
3365	}
3366	pluslist[jj-BO[7][1]]=BO[4][jj]+testid;
3367	}
3368	if(dontdrop[ii]==1)
3369	{
3370	ii++;
3371	continue;
3372	}
3373	if(!normalCross(pluslist))
3374	{
3375	//--- unfortunately, it is not normal crossing
3376	dontdrop[ii]=1;
3377	}
3378	}
3379	}
3380	//---------------------------------------------------------------------------
3381	// The returned list should look like the truncated output of inters_E
3382	//---------------------------------------------------------------------------
3383	list retlist;
3384	for(ii=1;ii<=size(E[4]);ii++)
3385	{
3386	if(dontdrop[ii]==1)
3387	{
3388	if(size(center)>0)
3389	{
3390	center=intersect(center,E[4][ii]);
3391	}
3392	else
3393	{
3394	center=E[4][ii];
3395	}
3396	}
3397	}
3398	retlist[1]=center;
3399	retlist[2]=E[2];
3400	retlist[3]=E[3];
3401	return(retlist);
3402	}
3403	//////////////////////////////////////////////////////////////////////////////
3404	static proc transversalT(ideal J, list E,list #)
3405	"Internal procedure - no help and no example available
3406	"
3407	{
3408	//----------------------------------------------------------------------------
3409	// check whether J and each element of the list E meet transversally
3410	//----------------------------------------------------------------------------
3411	def R=basering;
3412	if(size(#)>0)
3413	{
3414	ideal pp=#[1];
3415	}
3416	int i;
3417	ideal T,M;
3418	ideal Jstd=std(J);
3419	ideal Tstd;
3420	int d=nvars(basering)-dim(Jstd)+1; // d=n-dim(V(J) \cap hypersurface)
3421	for(i=1;i<=size(E);i++)
3422	{
3423	if(size(reduce(E[i],Jstd))==0)
3424	{
3425	//--- V(J) is contained in E[i]
3426	return(0);
3427	}
3428	T=J,E[i];
3429	Tstd=std(T);
3430	d=nvars(basering)-dim(Tstd);
3431	if(deg(Tstd[1])!=0)
3432	{
3433	//--- intersection is non-empty
3434	//!!! abgeklemmt, da es doch in der Praxis vorkommt und korrekt sein kann!!!
3435	//!!! wenn ueberhaupt dann -1 zurueckgeben!!!
3436	// if((d>=4)&&(size(T)>=10)){return(0);}
3437	M=minor(jacob(T),d,Tstd)+T;
3438	M=std(M);
3439	if(deg(M[1])>0)
3440	{
3441	//--- intersection is not transversal
3442	if(size(#)==0)
3443	{
3444	return(0);
3445	}
3446	M=std(radical(M));
3447	if(size(reduce(pp,M))>0){return(0);}
3448	}
3449	}
3450	}
3451	//--- passed all tests
3452	return(1);
3453	}
3454	///////////////////////////////////////////////////////////////////////////////
3455	static proc transversalTB(ideal J, list E,ideal V)
3456	"Internal procedure - no help and no example available
3457	"
3458	{
3459	//----------------------------------------------------------------------------
3460	// check whether J and each element of the list E meet transversally
3461	//----------------------------------------------------------------------------
3462	def R=basering;
3463
3464	int i;
3465	ideal T,M;
3466	ideal Jstd=std(J);
3467	ideal Tstd;
3468	int d=nvars(basering)-dim(Jstd)+1; // d=n-dim(V(J) \cap hypersurface)
3469	for(i=1;i<=size(E);i++)
3470	{
3471	if(size(reduce(E[i],Jstd))==0)
3472	{
3473	//--- V(J) is contained in E[i]
3474	return(0);
3475	}
3476	T=J,E[i];
3477	Tstd=std(T);
3478	d=nvars(basering)-dim(Tstd);
3479	if(deg(Tstd[1])!=0)
3480	{
3481	//--- intersection is non-empty
3482	if((d>=4)&&(size(T)>=10)){return(0);}
3483	M=minor(jacob(T),d,Tstd)+T;
3484	M=std(M+V);
3485	if(deg(M[1])>0)
3486	{
3487	return(0);
3488	}
3489	}
3490	}
3491	//--- passed all tests
3492	return(1);
3493	}
3494	///////////////////////////////////////////////////////////////////////////////
3495	static proc powerI(ideal I,int n,int m)
3496	{
3497	//--- compute (n!/m)-th power of I, more efficient variant
3498	int i;
3499	int mon=1;
3500	for(i=1;i<=size(I);i++)
3501	{
3502	if(size(I[i])>1){mon=0;break;}
3503	}
3504	if(mon)
3505	{
3506	if(size(reduce(I,std(radical(I[1]))))<size(I)-1){mon=0;}
3507	}
3508	if((mon)&&(size(I)>3))
3509	{
3510	int e=int(factorial(n))/m;
3511	ideal J=1;
3512	poly p=I[1];
3513	I=I[2..size(I)];
3514	ideal K=p^e;
3515	for(i=1;i<=e;i++)
3516	{
3517	J=interred(J*I);
3518	K=K,(p^(e-i))*J;
3519	}
3520	return(K);
3521	}
3522	for(i=n;i>1;i--)
3523	{
3524	if(i!=m)
3525	{
3526	I=I^i;
3527	}
3528	}
3529	return(I);
3530	}
3531
3532	///////////////////////////////////////////////////////////////////////////////
3533	static proc Coeff(list BO, int b, list #)
3534	"USAGE: Coeff (BO);
3535	@* BO = basic object, a list: ideal W,
3536	@* ideal J,
3537	@* intvec b (already truncated for Coeff),
3538	@* list Ex (already truncated for Coeff),
3539	@* ideal ab,
3540	@* intvec v,
3541	@* intvec w (already truncated for Coeff),
3542	@* matrix M
3543	@* b = integer indication bmax(BO)
3544	ASSUME: R = basering, a polynomial ring, W an ideal of R,
3545	@* J = ideal containing W
3546	COMPUTE: Coeff-Ideal of BO as defined in [Bravo,Encinas,Villamayor]
3547	RETURN: basic object of the Coeff-Ideal
3548	EXAMPLE: example Coeff; shows an example
3549	"
3550	{
3551	//!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
3552	//!!! TASK: lower dimension by more than one in a single step if possible !!!
3553	//!!! (improve bookkeeping of invariants in Coeff and Center) !!!
3554	//!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
3555	//---------------------------------------------------------------------------
3556	// Initialization and sanity checks
3557	//---------------------------------------------------------------------------
3558	int i,k,dummy,errtype;
3559	int ma=size(BO[4]);
3560	intvec merk;
3561	if(!defined(debugCoeff))
3562	{
3563	int debugCoeff;
3564	}
3565	ideal C;
3566	list L;
3567	if(size(#)!=0)
3568	{
3569	if(typeof(#[1])=="ideal")
3570	{
3571	L=#;
3572	}
3573	else
3574	{
3575	ma=#[1];
3576	L=DeltaList(BO);
3577	}
3578	}
3579	else
3580	{
3581	L=DeltaList(BO);
3582	}
3583
3584	if(debugCoeff)
3585	{
3586	"----> In Coeff: result of DeltaList:";
3587	L;
3588	}
3589	int d=size(L); // bmax of BO
3590	if((debugCoeff)&&(d!=b))
3591	{
3592	"!!!!!! Length of DeltaList does not equal second argument !!!!!!";
3593	"!!!!!! BO might not have been ord ~ 1 or wrong b !!!!!!";
3594	}
3595	if(b>=6){return(0);} // b is too big
3596	int e=int(factorial(b)); // b of Coeff-Ideal
3597	if(e==0)
3598	{
3599	ERROR( "// integer size too small for forming b! .");
3600	}
3601	if(b==0)
3602	{
3603	ERROR( "// second argument to Coeff should never be zero." );
3604	}
3605	//----------------------------------------------------------------------------
3606	// Form the Coeff-Ideal
3607	// Step 1: choose hypersurface
3608	// Step 2: sum over correct powers of Delta^i(BO[2])
3609	// Step 3: do the intersection
3610	//----------------------------------------------------------------------------
3611	//--- Step 1
3612	ideal Z;
3613	poly p;
3614	for(i=1;i<=ncols(L[d]);i++)
3615	{
3616	//--- Look for smooth hypersurface in generators of Delta^(bmax-1)(BO[2])
3617	dummy=goodChoice(BO,L[d][i]);
3618	if(!dummy)
3619	{
3620	Z= L[d][i];
3621	break;
3622	}
3623	else
3624	{
3625	if(dummy>1)
3626	{
3627	merk[size(merk)+1]=i;
3628	}
3629	if(dummy>errtype)
3630	{
3631	errtype=dummy;
3632	}
3633	}
3634	}
3635	if(size(Z)==0)
3636	{
3637	//--- no suitable element in generators of Delta^(bmax-1)(BO[2])
3638	//--- try random linear combination
3639	for(k=1;k<=10;k++)
3640	{
3641	for(i=2;i<=size(merk);i++)
3642	{
3643	p=p+random(-100,100)*L[d][merk[i]];
3644	}
3645	dummy=goodChoice(BO,p);
3646	if(!dummy)
3647	{
3648	Z=p;
3649	break;
3650	}
3651	else
3652	{
3653	p=0;
3654	}
3655	}
3656	if(dummy)
3657	{
3658	for(i=1;i<=size(L[d]);i++)
3659	{
3660	p=p+random(-100,100)*L[d][i];
3661	}
3662	dummy=goodChoice(BO,p);
3663	if(!dummy)
3664	{
3665	//--- found a suitable one
3666	Z=p;
3667	}
3668	}
3669	if(dummy)
3670	{
3671	//--- did not find a suitable random linear combination either
3672	if(dummy>errtype)
3673	{
3674	errtype=dummy;
3675	}
3676	list retlist=errtype,L[d];
3677	return(retlist);
3678	}
3679	}
3680	if(debugCoeff)
3681	{
3682	"----> In Coeff: Chosen hypersurface";
3683	Z;
3684	}
3685	//--- Step 2
3686	C=Z;
3687	for(i=0;i<b;i++)
3688	{
3689	C=C,powerI(simplify(reduce(L[i+1],std(Z)),2),b,b-i);
3690	}
3691	C=interred(C);
3692
3693	if(debugCoeff)
3694	{
3695	"----> In Coeff: J before saturation";
3696	C;
3697	}
3698
3699	//--- Step 3
3700	BO[1]=BO[1]+Z;
3701	BO[2]=C;
3702	for(i=1;i<=size(BO[4]);i++)
3703	{
3704	BO[6][i]=0; // reset intersection indicator
3705	BO[4][i]=BO[4][i]+Z; // intersect the E_i
3706	if(i<=ma)
3707	{
3708	BO[2]=sat(BO[2],BO[4][i]+BO[1])[1];
3709	// "strict transform" of J w.r.t E, not "total"
3710	}
3711	}
3712	if(debugCoeff)
3713	{
3714	"----> In Coeff:";
3715	" J after saturation:";
3716	BO[2];
3717	}
3718	return(BO);
3719	}
3720	example
3721	{"EXAMPLE:";
3722	echo = 2;
3723	ring R=0,(x,y,z),dp;
3724
3725	ideal W;
3726	ideal J=z^2+x^2*y^2;
3727	intvec b=0;
3728	list E;
3729	ideal abb=maxideal(1);
3730	intvec v;
3731	intvec w=-1;
3732	matrix M;
3733
3734	list BO=W,J,b,E,abb,v,w,M;
3735
3736	Coeff(BO,2);
3737	}
3738	//////////////////////////////////////////////////////////////////////////////
3739	static proc goodChoice(list BO, poly p)
3740	"Internal procedure - no help and no example available
3741	"
3742	{
3743	//---------------------------------------------------------------------------
3744	// test whether new W is smooth
3745	//---------------------------------------------------------------------------
3746	ideal W=BO[1]+ideal(p);
3747	if(size(reduce(p,std(BO[1])))==0)
3748	{
3749	//--- p is already in BO[1], i.e. does not define a hypersurface in W
3750	return(1);
3751	}
3752	if(dim(std(slocusE(W)))>=0)
3753	// if(dim(timeStd(slocusE(W),20))>=0)
3754	{
3755	//--- new W would not be smooth
3756	return(1);
3757	}
3758	if(size(BO[4])==0)
3759	{
3760	//--- E is empty, no further tests necessary
3761	return(0);
3762	}
3763	//--------------------------------------------------------------------------
3764	// test whether the hypersurface meets the E_i transversally
3765	//--------------------------------------------------------------------------
3766	list E=BO[4];
3767	int i,d;
3768	ideal T=W;
3769	ideal Tstd=std(T);
3770	d=nvars(basering)-dim(Tstd)+1;
3771	ideal M;
3772	for(i=1;i<=size(E);i++)
3773	{
3774	T=W,E[i];
3775	M=minor(jacob(T),d,Tstd)+T;
3776	M=std(M);
3777	if(deg(M[1])>0)
3778	{
3779	//--- intersection not transversal
3780	return(2);
3781	}
3782	}
3783	//--------------------------------------------------------------------------
3784	// test whether the new E_i have normal crossings
3785	//--------------------------------------------------------------------------
3786	for(i=1;i<=size(E);i++)
3787	{
3788	E[i]=E[i],p;
3789	}
3790	if(normalCross(E))
3791	{
3792	return(0);
3793	}
3794	else
3795	{
3796	return(2);
3797	}
3798	}
3799	//////////////////////////////////////////////////////////////////////////////
3800
3801	proc presentTree(list L)
3802	"USAGE: presentTree(L);
3803	L=list, output of resolve
3804	RETURN: nothing, only pretty printing of the output data of resolve()
3805	EXAMPLE: none
3806	"
3807	{
3808	def r=basering;
3809	int i,j,k;
3810	if(size(L[2])==1)
3811	{
3812	"The object was already resolved or the list L does not";
3813	"have required input format. There is just one chart in";
3814	"the tree.";
3815	return();
3816	}
3817	for(i=1;i<=size(L[1]);i++)
3818	{
3819	" ";
3820	"/////////////////////////// Final Chart",i,"/////////////////////////";
3821	def s=L[1][i];
3822	setring s;
3823	"======================== History of this chart ======================";
3824	for(j=2;j<=ncols(path);j++)
3825	{
3826	" ";
3827	"Blow Up",j-1,":";
3828	" Center determined in L[2]["+string(path[1,j])+"],";
3829	" Passing to chart ",path[2,j]," in resulting blow up.";
3830	}
3831	" ";
3832	"======================== Data of this chart ========================";
3833	showBO(BO);
3834	setring r;
3835	kill s;
3836	pause();
3837	}
3838	"////////////////////////////////////////////////////////////////////";
3839	"For identification of exceptional divisors please use the tools";
3840	"provided by reszeta.lib, e.g. collectDiv.";
3841	"For viewing an illustration of the tree of charts please use the";
3842	"procedure ResTree from resgraph.lib.";
3843	"////////////////////////////////////////////////////////////////////";
3844	return();
3845	}
3846	//////////////////////////////////////////////////////////////////////////////
3847
3848	proc showBO(list BO)
3849	"USAGE: showBO(BO);
3850	@* BO=basic object, a list: ideal W,
3851	@* ideal J,
3852	@* intvec b (already truncated for Coeff),
3853	@* list Ex (already truncated for Coeff),
3854	@* ideal ab,
3855	@* intvec v,
3856	@* intvec w (already truncated for Coeff),
3857	@* matrix M
3858	RETURN: nothing, only pretty printing
3859	EXAMPLE: none
3860	"
3861	{
3862	" ";
3863	"==== Ambient Space: ";BO[1];" ";
3864	"==== Ideal of Variety: ";BO[2];" ";
3865	int i;
3866	list M;
3867	for(i=1;i<=size(BO[4]);i++)
3868	{
3869	M[i]=ideal(BO[4][i]);
3870	}
3871	"==== Exceptional Divisors: ";print(M);" ";
3872	"==== Images of variables of original ring:";BO[5];" ";
3873	}
3874	//////////////////////////////////////////////////////////////////////////////
3875	//////////////////////// main procedure ////////////////////////////////
3876	//////////////////////////////////////////////////////////////////////////////
3877	proc resolve(ideal J, list #)
3878	"USAGE: resolve (J); or resolve (J,i[,k]);
3879	@* J ideal
3880	@* i,k int
3881	COMPUTE: a resolution of J,
3882	@* if i > 0 debugging is turned on according to the following switches:
3883	@* j1: value 0 or 1; turn off or on correctness checks in all steps
3884	@* j2: value 0 or 2; turn off or on debugCenter
3885	@* j3: value 0 or 4; turn off or on debugBlowUp
3886	@* j4: value 0 or 8; turn off or on debugCoeff
3887	@* j5: value 0 or 16:turn off or on debugging of Intersection with E^-
3888	@* j6: value 0 or 32:turn off or on stop after pass throught the loop
3889	@* i=j1+j2+j3+j4+j5+j6
3890	RETURN: a list l of 2 lists of rings
3891	l[1][i] is a ring containing a basic object BO, the result of the
3892	resolution.
3893	l[2] contains all rings which occured during the resolution process
3894	NOTE: result may be viewed in a human readable form using presentTree()
3895	EXAMPLE: example resolve; shows an example
3896	"
3897	{
3898	//----------------------------------------------------------------------------
3899	// Initialization and sanity checks
3900	//----------------------------------------------------------------------------
3901	def R=basering;
3902	list allRings;
3903	allRings[1]=R;
3904	list endRings;
3905	module path=[0,-1];
3906	ideal W;
3907	list E;
3908	ideal abb=maxideal(1);
3909	intvec v;
3910	intvec bvec;
3911	intvec w=-1;
3912	matrix intE;
3913	int extra,bm;
3914	if(defined(BO)){kill BO;}
3915	if(defined(cent)){kill cent;}
3916
3917	ideal Jrad=equiRadical(J);
3918	if(size(reduce(Jrad,std(J)))!=0)
3919	{
3920	"WARNING! The input is not reduced or not equidimensional!";
3921	"We will continue with the reduced top-dimensional part of input";
3922	J=Jrad;
3923	}
3924
3925	int i,j,debu,loca,locaT,ftemp,debugResolve,smooth;
3926	//--- switches for local and for debugging may occur in any order
3927	i=size(#);
3928	extra=3;
3929	for(j=1;j<=i;j++)
3930	{
3931	if(typeof(#[j])=="int")
3932	{
3933	debugResolve=#[j];
3934	//--- debu: debug switch for resolve, smallest bit in debugResolve
3935	debu=debugResolve mod 2;
3936	}
3937	else
3938	{
3939	if(#[j]=="M")
3940	{
3941	bm=1;
3942	ERROR("Not implemented yet");
3943	}
3944	if(#[j]=="E"){extra=0;}
3945	if(#[j]=="A"){extra=2;}
3946	if(#[j]=="K"){extra=3;}
3947	if(#[j]=="L"){loca=1;}
3948	}
3949	}
3950	if(loca)
3951	{
3952	list qs=minAssGTZ(J);
3953	ideal K=ideal(1);
3954	for(j=1;j<=size(qs);j++)
3955	{
3956	if(size(reduce(qs[j],std(maxideal(1))))==0)
3957	{
3958	K=intersect(K,qs[j]);
3959	}
3960	}
3961	J=K;
3962	list qr=minAssGTZ(slocus(J));
3963	K=ideal(1);
3964	for(j=1;j<=size(qr);j++)
3965	{
3966	if(size(reduce(qr[j],std(maxideal(1))))!=0)
3967	{
3968	K=intersect(K,qr[j]);
3969	smooth++;
3970	}
3971	else
3972	{
3973	if(dim(std(qr[j]))>0){loca=0;}
3974	//---- test for isolated singularity at 0
3975	}
3976	}
3977	K=std(K);
3978	//---- if deg(K[1])==0 the point 0 is on all components of the singular
3979	//---- locus and we can work globally
3980	if(smooth==size(qr)){smooth=-1;}
3981	//---- the point 0 is not on the singular locus
3982	if((deg(K[1])>0)&&(smooth>=0)&&(!loca))
3983	{
3984	locaT=1;
3985	poly @p;
3986	for(j=1;j<=size(K);j++)
3987	{
3988	if(jet(K[j],0)!=0)
3989	{
3990	@p=K[j];
3991	break;
3992	}
3993	}
3994	export(@p);
3995	}
3996	if((loca)&&(!smooth)){loca=0;}
3997	//---- the case that 0 is isolated singularity and the only singular point
3998	}
3999	export(locaT);
4000	//---In case of option "L" the following holds
4001	//---loca=0 and locaT=0 we perform the global case
4002	//---loca !=0: 0 is isolated singular point, but there are other singularities
4003	//---locaT!=0: O is singular point, but not isolated, and there is a componente//--- of the singular locus not containing 0
4004
4005	//--- if necessary, set the corresponding debugFlags
4006	if(defined(debugResolve))
4007	{
4008	//--- 2nd bit from the right
4009	int debugCenter=(debugResolve div 2) mod 2;
4010	export debugCenter;
4011	//--- 3rd bit from the right
4012	int debugBlowUp=(debugResolve div 4) mod 2;
4013	export debugBlowUp;
4014	//--- 4th bit from the right
4015	int debugCoeff=(debugResolve div 8) mod 2;
4016	export debugCoeff;
4017	//--- 5th bit from the right
4018	int debug_Inters_E=(debugResolve div 16) mod 2;
4019	export debug_Inters_E;
4020	//--- 6th bit from the right
4021	int praes_stop=(debugResolve div 32) mod 2;
4022	}
4023	//--- set the correct attributes to J for speed ups
4024	if( typeof(attrib(J,"isEqui"))!="int" )
4025	{
4026	if(size(J)==1)
4027	{
4028	attrib(J,"isEqui",1);
4029	}
4030	else
4031	{
4032	attrib(J,"isEqui",0);
4033	}
4034	}
4035	if(size(J)==1)
4036	{
4037	attrib(J,"isHy",1);
4038	}
4039	else
4040	{
4041	attrib(J,"isHy",0);
4042	}
4043	//--- create the BO
4044	list BO=W,J,bvec,E,abb,v,w,intE;
4045	if(defined(invSat)){kill invSat;}
4046	list invSat=ideal(0),intvec(0);
4047	export(invSat);
4048	if(bm)
4049	{
4050	intmat invmat[2][1]=0,-1;
4051	BO[9]=invmat;
4052	}
4053	else
4054	{
4055	BO[9]=intvec(0);
4056	}
4057	export BO;
4058	list tmpList;
4059	int blo;
4060	int k,Ecount,tmpPtr;
4061	i=0;
4062	if(smooth==-1)
4063	{
4064	endRings[1]=R;
4065	list result=endRings,allRings;
4066	if(debu)
4067	{
4068	"============= result will be tested ==========";
4069	" ";
4070	"the number of charts obtained:",size(endRings);
4071	"============= result is o.k. ==========";
4072	}
4073	kill debugCenter,debugBlowUp,debugCoeff,debug_Inters_E;
4074	return(result);
4075	}
4076	//-----------------------------------------------------------------------------
4077	// While there are rings to be considered, determine center and blow up
4078	//-----------------------------------------------------------------------------
4079	while(i<size(allRings))
4080	{
4081	i++;
4082	def S=allRings[i];
4083	setring S;
4084	list pr;
4085	ideal Jstd=std(BO[2]);
4086	//-----------------------------------------------------------------------------
4087	// Determine Center
4088	//-----------------------------------------------------------------------------
4089	if(i==1)
4090	{
4091	list deltaL=DeltaList(BO);
4092	ideal sL=radical(deltaL[size(deltaL)]);
4093	if((deg(std(slocus(sL))[1])==0)&&(size(minAssGTZ(sL))==1))
4094	{
4095	list @ce=sL,intvec(-1),intvec(0),intvec(0);
4096	ideal cent=@ce[1];
4097	}
4098	}
4099	//--- before computing a center, check whether we have a stored one
4100	if(size(BO)>9)
4101	{
4102	while(size(BO[10])>0)
4103	{
4104	list @ce=BO[10][1];
4105	//--- check of the center
4106	// @ce=correctC(BO,@ce,bm);
4107	//--- use stored center
4108	BO[10]=delete(BO[10],1);
4109	if(size(@ce[1])==0)
4110	{
4111	//--- stored center was not ok
4112	continue;
4113	}
4114	tmpPtr=0;
4115	for(Ecount=1;Ecount <= size(@ce[2]); Ecount++)
4116	{
4117	if(@ce[2][Ecount]>-1)
4118	{
4119	tmpPtr=tmpPtr+@ce[2][Ecount];
4120	}
4121	else
4122	{
4123	@ce[2][Ecount]=size(BO[4])-tmpPtr-1;
4124	for(int cnthlp=1;cnthlp<=size(BO[10]);cnthlp++)
4125	{
4126	BO[10][cnthlp][2][Ecount]=@ce[2][Ecount];
4127	}
4128	kill cnthlp;
4129	break;
4130	}
4131	}
4132	if(Ecount<size(@ce[2]))
4133	{
4134	for(tmpPtr=Ecount+1;tmpPtr<=size(@ce[2]);tmpPtr++)
4135	{
4136	@ce[2][tmpPtr]=0;
4137	for(int cnthlp=1;cnthlp<=size(BO[10]);cnthlp++)
4138	{
4139	BO[10][cnthlp][2][tmpPtr]=@ce[2][tmpPtr];
4140	}
4141	kill cnthlp;
4142	}
4143	}
4144	break;
4145	}
4146	if(defined(@ce))
4147	{
4148	if(size(@ce[1])==0)
4149	{
4150	kill @ce;
4151	}
4152	else
4153	{
4154	ideal cent=@ce[1];
4155	}
4156	}
4157	if(size(BO[10])==0)
4158	{
4159	//--- previously had stored centers, all have been used; we need to clean up
4160	BO=delete(BO,10);
4161	}
4162	}
4163	if((loca)&&(i==1))
4164	{
4165	//--- local case: initial step is blow-up in origin
4166	if(defined(@ce)){kill @ce;}
4167	if(defined(cent)){kill cent;}
4168	if(size(reduce(slocusE(BO[2]),std(maxideal(1))))==0)
4169	{
4170	list @ce=maxideal(1),intvec(-1),intvec(0),intvec(0);
4171	}
4172	else
4173	{
4174	list @ce=BO[2],intvec(-1),intvec(1),intvec(0);
4175	}
4176	ideal cent=@ce[1];
4177	}
4178	if(((loca)\|\|(locaT))&&(i!=1))
4179	{
4180	int JmeetsE;
4181	for(j=1;j<=size(BO[4]);j++)
4182	{
4183	if(deg(std(BO[2]+BO[4][j])[1])!=0)
4184	{
4185	JmeetsE=1;
4186	break;
4187	}
4188	}
4189	if(!JmeetsE)
4190	{
4191	list @ce=BO[2],intvec(-1),intvec(1),intvec(0);
4192	ideal cent=@ce[1];
4193	}
4194	kill JmeetsE;
4195	}
4196	if((locaT)&&(!defined(@ce)))
4197	{
4198	if(@p!=1)
4199	{
4200	list tr=minAssGTZ(slocusE(BO[2]));
4201	ideal L=ideal(1);
4202	for(j=1;j<=size(tr);j++)
4203	{
4204	if(size(reduce(ideal(@p),std(tr[j])))==0)
4205	{
4206	L=intersect(L,tr[j]);
4207	}
4208	}
4209	L=std(L);
4210	if(deg(L[1])==0)
4211	{
4212	@p=1;
4213	}
4214	else
4215	{
4216	ideal fac=factorize(@p,1);
4217	if(size(fac)==1)
4218	{
4219	@p=fac[1];
4220	}
4221	else
4222	{
4223	for(j=1;j<=size(fac);j++)
4224	{
4225	if(reduce(fac[j],L)==0)
4226	{
4227	@p=fac[j];
4228	break;
4229	}
4230	}
4231	}
4232	}
4233	kill tr,L;
4234	}
4235	execute("ring R1=("+charstr(S)+"),(@z,"+varstr(S)+"),dp;");
4236	poly p=imap(S,@p);
4237	list BO=imap(S,BO);
4238	list invSat=imap(S,invSat);
4239	export(invSat);
4240	ideal te=std(BO[2]);
4241	BO[1]=BO[1]+ideal(@z*p-1);
4242	BO[2]=BO[2]+ideal(@z*p-1);
4243	for(j=1;j<=size(BO[4]);j++)
4244	{
4245	BO[4][j]=BO[4][j]+ideal(@z*p-1);
4246	}
4247	//--- for computation of center: drop components not meeting the Ei
4248	def BO2=BO;
4249	list qs=minAssGTZ(BO2[2]);
4250	ideal K=ideal(1);
4251	for(j=1;j<=size(qs);j++)
4252	{
4253	if(CompMeetsE(qs[j],BO2[4]))
4254	{
4255	K=intersect(K,qs[j]);
4256	}
4257	}
4258	BO2[2]=K;
4259	//--- check whether we are done
4260	if(deg(std(BO2[2])[1])==0)
4261	{
4262	list @ce=BO[2],intvec(-1),intvec(1),intvec(0);
4263	}
4264	if(!defined(@ce))
4265	{
4266	if(bm)
4267	{
4268	list @ce=CenterBM(BO2);
4269	}
4270	else
4271	{
4272	list @ce=CenterBO(BO2);
4273	}
4274	}
4275	//--- if computation of center returned BO2[2], we are done
4276	//--- ==> set @ce to BO[2], because later checks work with BO instead of BO2
4277	if((size(reduce(@ce[1],std(BO2[2])))==0)&&
4278	(size(reduce(BO2[2],std(@ce[1])))==0))
4279	{
4280	@ce[1]=BO[2];
4281	}
4282	if(size(specialReduce(@ce[1],te,p))==0)
4283	{
4284	BO=imap(S,BO);
4285	@ce[1]=BO[2];
4286	}
4287	else
4288	{
4289	//@ce=correctC(BO,@ce,bm);
4290	@ce[1]=eliminate(@ce[1],@z);
4291	}
4292	setring S;
4293	list @ce=imap(R1,@ce);
4294	kill R1;
4295
4296	if((size(reduce(BO[2],std(@ce[1])))==0)
4297	&&(size(reduce(@ce[1],Jstd))==0))
4298	{
4299	//--- J and center coincide
4300	pr[1]=@ce[1];
4301	ideal cent=@ce[1];
4302	}
4303	else
4304	{
4305	//--- decompose center and use first component
4306	pr=minAssGTZ(@ce[1]);
4307	if(size(reduce(@p,std(pr[1])))==0){"Achtung";~;}
4308	if(deg(std(slocus(pr[1]))[1])>0){"singulaer";~;}
4309	ideal cent=pr[1];
4310	}
4311	if(size(pr)>1)
4312	{
4313	//--- store the other components
4314	for(k=2;k<=size(pr);k++)
4315	{
4316	if(size(reduce(@p,std(pr[k])))==0){"Achtung";~;}
4317	if(deg(std(slocus(pr[k]))[1])>0){"singulaer";~;}
4318	if(size(reduce(@p,std(pr[k])))!=0)
4319	{
4320	tmpList[size(tmpList)+1]=list(pr[k],@ce[2],@ce[3],@ce[4]);
4321	}
4322	}
4323	BO[10]=tmpList;
4324	kill tmpList;
4325	list tmpList;
4326	}
4327	}
4328	if(!defined(@ce))
4329	{
4330	//--- no previously determined center, we need to compute one
4331	if(loca)
4332	{
4333	//--- local case: center should be inside exceptional locus
4334	ideal Ex=ideal(1);
4335	k=0;
4336	for(j=1;j<=size(BO[4]);j++)
4337	{
4338	if(deg(BO[4][j][1])!=0)
4339	{
4340	Ex=Ex*BO[4][j]; //----!!!!hier evtl. Durchschnitt???
4341	k++;
4342	}
4343	}
4344	//--- for computation of center: drop components not meeting the Ei
4345	list BOloc=BO;
4346	list qs=minAssGTZ(BOloc[2]);
4347	ideal K=ideal(1);
4348	for(j=1;j<=size(qs);j++)
4349	{
4350	if(CompMeetsE(qs[j],BOloc[4]))
4351	{
4352	K=intersect(K,qs[j]);
4353	}
4354	}
4355	BOloc[2]=K;
4356	//--- check whether we are done
4357	if(deg(std(BOloc[2])[1])==0)
4358	{
4359	list @ce=BO[2],intvec(-1),intvec(1),intvec(0);
4360	}
4361	if(!defined(@ce))
4362	{
4363	if(BO[3][1]!=0)
4364	{
4365	BOloc[2]=BO[2]+Ex^((BO[3][1] div k)+1);//!!!!Vereinfachen???
4366	}
4367	else
4368	{
4369	BOloc[2]=BO[2]+Ex^((size(DeltaList(BO)) div k)+1);
4370	}
4371	if(bm)
4372	{
4373	list @ce=CenterBM(BOloc);
4374	}
4375	else
4376	{
4377	list @ce=CenterBO(BOloc);
4378	}
4379	if(size(reduce(Ex,std(@ce[1])))!=0)
4380	{
4381	list tempPr=minAssGTZ(@ce[1]);
4382	for(k=size(tempPr);k>=1;k--)
4383	{
4384	if(size(reduce(Ex,std(tempPr[k])))!=0)
4385	{
4386	tempPr=delete(tempPr,k);
4387	}
4388	}
4389	@ce[1]=1;
4390	for(k=1;k<=size(tempPr);k++)
4391	{
4392	@ce[1]=intersect(@ce[1],tempPr[k]);
4393	}
4394	if(deg(std(@ce[1])[1])==0)
4395	{
4396	@ce[1]=BO[2];
4397	}
4398	}
4399	}
4400	//--- test whether we are done
4401	if(size(reduce(slocusE(BO[2]),std(@ce[1])))!=0)
4402	{
4403	if(transversalT(BO[2],BO[4]))
4404	{
4405	if(defined(E)){kill E;}
4406	list E=BO[4];
4407	for(j=1;j<=size(E);j++){if(deg(E[j][1])>0){E[j]=E[j]+BO[2];}}
4408	if(normalCross(E))
4409	{
4410	@ce[1]=BO[2];
4411	}
4412	kill E;
4413	}
4414	}
4415	}
4416	else
4417	{
4418	//--- non-local
4419	if(bm)
4420	{
4421	list @ce=CenterBM(BO);
4422	}
4423	else
4424	{
4425	list @ce=CenterBO(BO);
4426	}
4427	//--- check of the center
4428	//@ce=correctC(BO,@ce,bm);
4429	if((size(@ce[1])==0)&&(size(@ce[4])<(size(@ce[3])-1)))
4430	{
4431	intvec xxx=@ce[3];
4432	xxx=xxx[1..size(@ce[4])];
4433	@ce[3]=xxx;
4434	xxx=@ce[2];
4435	xxx=xxx[1..size(@ce[4])];
4436	@ce[2]=xxx;
4437	kill xxx;
4438	}
4439	}
4440	if((size(reduce(BO[2],std(@ce[1])))==0)
4441	&&(size(reduce(@ce[1],Jstd))==0))
4442	{
4443	//--- J and center coincide
4444	pr[1]=@ce[1];
4445	ideal cent=@ce[1];
4446	}
4447	else
4448	{
4449	//--- decompose center and use first component
4450	pr=minAssGTZ(@ce[1]);
4451	ideal cent=pr[1];
4452	}
4453	if(size(pr)>1)
4454	{
4455	//--- store the other components
4456	for(k=2;k<=size(pr);k++)
4457	{
4458	tmpList[k-1]=list(pr[k],@ce[2],@ce[3],@ce[4]);
4459	}
4460	BO[10]=tmpList;
4461	kill tmpList;
4462	list tmpList;
4463	}
4464	}
4465	//--- do not forget to update BO[7] and BO[3]
4466	export cent;
4467	BO[7]=@ce[2];
4468	BO[3]=@ce[3];
4469	if((loca\|\|locaT)&&(size(@ce)<4)){@ce[4]=0;} //Provisorium !!!
4470	if((size(@ce[4])<size(@ce[2])-1)\|\|(size(@ce[4])<size(@ce[3])-1))
4471	{
4472	if((deg(std(@ce[1])[1])==0)&&(deg(std(BO[2])[1])==0))
4473	{
4474	intvec nullvec;
4475	nullvec[size(@ce[2])-1]=0;
4476	@ce[4]=nullvec;
4477	kill nullvec;
4478	}
4479	else
4480	{
4481	"ERROR:@ce[4] hat falsche Laenge - nicht-trivialer Fall";
4482	~;
4483	}
4484	}
4485	if((typeof(@ce[4])=="intvec") \|\| (typeof(@ce[4])=="intmat"))
4486	{
4487	BO[9]=@ce[4];
4488	}
4489	//---------------------------------------------------------------------------
4490	// various checks and debug output
4491	//---------------------------------------------------------------------------
4492	if((debu) \|\| (praes_stop))
4493	{
4494	//--- Show BO of this step
4495	"++++++++++++++ Overview of Current Chart +++++++++++++++++++++++";
4496	"Current number of final charts:",size(endRings);
4497	"Total number of charts currently in chart-tree:",size(allRings);
4498	"Index of the current chart in chart-tree:",i;
4499	"++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++";
4500	showBO(BO);
4501	"-------------------------- Upcoming Center ---------------------";
4502	interred(cent);
4503	"----------------------------------------------------------------";
4504	}
4505	if(praes_stop)
4506	{
4507	~;
4508	}
4509	if(debu)
4510	{
4511	//--- various checks, see output for comments
4512	if(size(BO[1])>0)
4513	{
4514	if(deg(BO[1][1])==0)
4515	{
4516	"!!! W is empty !!!";
4517	path;
4518	setring R;
4519	kill S;
4520	list result=endRings,allRings;
4521	return(result);
4522	}
4523	if(deg(std(slocusE(BO[1]))[1])>0)
4524	{
4525	"!!! W not smooth !!!";
4526	path;
4527	setring R;
4528	kill S;
4529	list result=endRings,allRings;
4530	return(result);
4531	}
4532	}
4533	if((!loca)&&(!locaT))
4534	{
4535	if(deg(std(slocusE(cent+BO[1]))[1])>0)
4536	{
4537	"!!! Center not smooth !!!";
4538	path;
4539	std(cent+BO[1]);
4540	~;
4541	setring R;
4542	kill S;
4543	list result=endRings,allRings;
4544	return(result);
4545	}
4546	}
4547	for(j=1;j<=size(BO[4]);j++)
4548	{
4549	if(deg(BO[4][j][1])>0)
4550	{
4551	if(deg(std(slocusE(BO[4][j]+BO[1]))[1])>0)
4552	{
4553	"!!! exceptional divisor is not smooth !!!";
4554	path;
4555	setring R;
4556	kill S;
4557	list result=endRings,allRings;
4558	return(result);
4559	}
4560	}
4561	}
4562	if((!loca)&&(!locaT))
4563	{
4564	if((norC(BO,cent))&&(size(reduce(cent,Jstd))!=0))
4565	{
4566	"!!! this chart is already finished !!!";
4567	cent=BO[2];
4568	~;
4569	}
4570	}
4571	}
4572	//----------------------------------------------------------------------------
4573	// Do the blow up
4574	//----------------------------------------------------------------------------
4575	//!!!! Change this as soon as there is time!!!
4576	//!!!! quick and dirty bug fix for old shortcut which has not yet been killed
4577	if((dim(std(cent))==0)&&defined(shortcut)) {kill shortcut;}
4578	//!!! end of bugfix
4579	if(size(reduce(cent,Jstd))!=0)
4580	{
4581	//--- center does not equal J
4582	tmpList=blowUpBO(BO,cent,extra);
4583	if((debu)&&(!loca)&&(!locaT))
4584	{
4585	//--- test it, if debu is set
4586	if(!testBlowUp(BO,cent,tmpList,i,extra))
4587	{
4588	"non-redundant chart has been killed!";
4589	~;
4590	}
4591	}
4592	//--- extend the list of all rings
4593	allRings[size(allRings)+1..size(allRings)+size(tmpList)]=
4594	tmpList[1..size(tmpList)];
4595	for(j=1;j<=size(tmpList);j++)
4596	{
4597	def Q=allRings[size(allRings)-j+1];
4598	setring Q;
4599	def path=imap(S,path);
4600	path=path,[i,size(tmpList)-j+1];
4601	export path;
4602	setring S;
4603	kill Q;
4604	}
4605	kill tmpList;
4606	list tmpList;
4607	}
4608	else
4609	{
4610	//--- center equals J
4611	k=0;
4612	for(j=1;j<=size(BO[6]);j++)
4613	{
4614	if(BO[6][j]!=1)
4615	{
4616	//--- there is an E_i which meets J in this chart
4617	k=1;
4618	break;
4619	}
4620	}
4621	if((k)\|\|(extra==2))
4622	{
4623	//--- chart finished, non-redundant
4624	endRings[size(endRings)+1]=S;
4625	}
4626	}
4627	kill pr;
4628	setring R;
4629	kill S;
4630	}
4631	//---------------------------------------------------------------------------
4632	// set up the result, test it (if debu is set) and return it
4633	//---------------------------------------------------------------------------
4634	list result=endRings,allRings;
4635	if(debu)
4636	{
4637	"============= result will be tested ==========";
4638	" ";
4639	"the number of charts obtained:",size(endRings);
4640	if(locaT){loca=2;}
4641	int tes=testRes(J,endRings,loca);
4642	if(tes)
4643	{
4644	"============= result is o.k. ==========";
4645	}
4646	else
4647	{
4648	"============ result is wrong =========="; ~;
4649	}
4650	}
4651	kill debugCenter,debugBlowUp,debugCoeff,debug_Inters_E;
4652	if(locaT){kill @p;}
4653	kill locaT;
4654	return(result);
4655	}
4656	example
4657	{"EXAMPLE:";
4658	echo = 2;
4659	ring R=0,(x,y,z),dp;
4660	ideal J=x3+y5+yz2+xy4;
4661	list L=resolve(J,0);
4662	def Q=L[1][7];
4663	setring Q;
4664	showBO(BO);
4665	}
4666	//////////////////////////////////////////////////////////////////////////
4667	//static
4668	proc CompMeetsE(ideal J, list E)
4669	"Internal procedure - no help and no example available
4670	"
4671	{
4672	int i;
4673	for(i=1;i<=size(E);i++)
4674	{
4675	if(deg(std(E[i])[1])!=0)
4676	{
4677	if(deg(std(J+E[i])[1])!=0)
4678	{
4679	return(1);
4680	}
4681	}
4682	}
4683	return(0);
4684	}
4685
4686	//========================================================================
4687	//-------------- procedures for testing the result ----------------------
4688	// (not yet commented)
4689	//========================================================================
4690
4691	//////////////////////////////////////////////////////////////////////////
4692	static proc testRes(ideal J,list L,int loca)
4693	"Internal procedure - no help and no example available
4694	"
4695	{
4696	int loc;
4697	if(defined(locaT)){loc=locaT;}
4698	if(loc){loca=0;}
4699	def R=basering;
4700	ideal M=maxideal(1);
4701	int i,j,tr;
4702	for(i=1;i<=size(L);i++)
4703	{
4704	def Q=L[i];
4705	setring Q;
4706	ideal J=BO[2];
4707	list E=BO[4];
4708	map phi=R,BO[5];
4709	ideal K=phi(J)+BO[1];
4710	ideal stTK=std(K);
4711
4712	if(loca)
4713	{
4714	ideal M=phi(M)+BO[1];
4715	ideal stTM=std(M);
4716	}
4717	for(j=1;j<=size(E);j++)
4718	{
4719	if(deg(E[j][1])>0)
4720	{
4721	stTK=sat(stTK,E[j])[1];
4722	}
4723	if(loca)
4724	{
4725	stTM=sat(stTM,E[j])[1];
4726	}
4727	}
4728	ideal sL=slocusE(J);
4729	if(loca){sL=sL+stTM;}
4730	ideal sLstd=std(sL);
4731	if(deg(sLstd[1])>0)
4732	{
4733	if(!loc)
4734	{
4735	"J is not smooth";i;
4736	setring R;
4737	return(0);
4738	}
4739	if(size(reduce(@p,std(radical(sLstd))))>0)
4740	{
4741	"J is not smooth";i;
4742	setring R;
4743	return(0);
4744	}
4745	}
4746	if(!((size(reduce(J,std(stTK)))==0)
4747	&&(size(reduce(stTK,std(J)))==0)))
4748	{
4749	"map is wrong";i;
4750	setring R;
4751	return(0);
4752	}
4753	if(loc){tr=transversalT(J,E,@p);}
4754	else{tr=transversalT(J,E);}
4755	if(!tr)
4756	{
4757	"E not transversal with J";i;
4758	setring R;
4759	return(0);
4760	}
4761	if(!normalCross(E))
4762	{
4763	"E not normal crossings";i;
4764	setring R;
4765	return(0);
4766	}
4767	for(j=1;j<=size(E);j++)
4768	{
4769	if(deg(E[j][1])>0){E[j]=E[j]+J;}
4770	}
4771	if(!normalCross(E))
4772	{
4773	"E not normal crossings with J";i;
4774	setring R;
4775	return(0);
4776	}
4777	kill J,E,phi,K,stTK;
4778	if(loca){kill M,stTM;}
4779	setring R;
4780	kill Q;
4781	}
4782	return(1);
4783	}
4784	//////////////////////////////////////////////////////////////////////////////
4785	static proc testBlowUp(list BO,ideal cent,list tmpList, int j, int extra)
4786	{
4787	def R=basering;
4788	int n=nvars(basering);
4789	int i;
4790	if((extra!=3)&&(extra!=2))
4791	{
4792	ideal K=BO[1],BO[2],cent;
4793	for(i=1;i<=size(BO[4]);i++)
4794	{
4795	K=K,BO[4][i];
4796	}
4797	list N=findvars(K);
4798	//list N=findvars(BO[2]);
4799	if(size(N[1])<n)
4800	{
4801	string newvar=string(N[1]);
4802	execute("ring R1=("+charstr(R)+"),("+newvar+"),dp;");
4803	list BO=imap(R,BO);
4804	ideal cent=imap(R,cent);
4805	n=nvars(R1);
4806	}
4807	else
4808	{
4809	def R1=basering;
4810	}
4811	}
4812	else
4813	{
4814	def R1=basering;
4815	}
4816
4817	i=0;
4818	ideal T=cent;
4819	ideal TW;
4820	for(i=1;i<=size(tmpList);i++)
4821	{
4822	def Q=tmpList[i];
4823	setring Q;
4824	map phi=R1,lastMap;
4825	ideal TE=radical(slocusE(BO[2]));
4826	setring R1;
4827	TW=preimage(Q,phi,TE);
4828	T=intersect(T,TW);
4829	kill Q;
4830	}
4831	ideal sL=intersect(slocusE(BO[2]),cent);
4832	if(size(reduce(sL,std(radical(T))))>0){setring R;return(0);}
4833	if(size(reduce(T,std(radical(sL))))>0){setring R;return(0);}
4834	setring R;
4835	return(1);
4836	}
4837	//////////////////////////////////////////////////////////////////////////////
4838	static proc normalCross(list E,list #)
4839	"Internal procedure - no help and no example available
4840	"
4841	{
4842	int loc;
4843	if((defined(locaT))&&(defined(@p)))
4844	{
4845	loc=1;
4846	ideal pp=@p;
4847	}
4848	int i,d,j;
4849	int n=nvars(basering);
4850	list E1,E2;
4851	ideal K,M,Estd,cent;
4852	intvec v,w;
4853	if(size(#)>0){cent=#[1];}
4854
4855	for(i=1;i<=size(E);i++)
4856	{
4857	Estd=std(E[i]);
4858	if(deg(Estd[1])>0)
4859	{
4860	E1[size(E1)+1]=Estd;
4861	}
4862	}
4863	E=E1;
4864	for(i=1;i<=size(E);i++)
4865	{
4866	v=i;
4867	E1[i]=list(E[i],v);
4868	}
4869	list ll;
4870	int re=1;
4871	int ok;
4872	while(size(E1)>0)
4873	{
4874	K=E1[1][1];
4875	v=E1[1][2];
4876	attrib(K,"isSB",1);
4877	E1=delete(E1,1);
4878	d=n-dim(K);
4879	M=minor(jacob(K),d)+K;
4880	if(deg(std(M)[1])>0)
4881	{
4882	if(size(#)>0)
4883	{
4884	if(size(reduce(M,std(cent)))>0)
4885	{
4886	ll[size(ll)+1]=std(M);
4887	}
4888	else
4889	{
4890	ok=1;
4891	}
4892	}
4893	if(!loc)
4894	{
4895	re=0;
4896	}
4897	else
4898	{
4899	if(size(reduce(pp,std(radical(M))))>0){re=0;}
4900	}
4901	}
4902	for(i=1;i<=size(E);i++)
4903	{
4904	for(j=1;j<=size(v);j++){if(v[j]==i){break;}}
4905	if(j<=size(v)){if(v[j]==i){i++;continue;}}
4906	Estd=std(K+E[i]);
4907	w=v;
4908	if(deg(Estd[1])==0){i++;continue;}
4909	if(d==n-dim(Estd))
4910	{
4911	if(size(#)>0)
4912	{
4913	if(size(reduce(Estd,std(cent)))>0)
4914	{
4915	ll[size(ll)+1]=Estd;
4916	}
4917	else
4918	{
4919	ok=1;
4920	}
4921	}
4922	if(!loc)
4923	{
4924	re=0;
4925	}
4926	else
4927	{
4928	if(size(reduce(pp,std(radical(M))))>0){re=0;}
4929	}
4930	}
4931	w[size(w)+1]=i;
4932	E2[size(E2)+1]=list(Estd,w);
4933	}
4934	if(size(E2)>0)
4935	{
4936	if(size(E1)>0)
4937	{
4938	E1[size(E1)+1..size(E1)+size(E2)]=E2[1..size(E2)];
4939	}
4940	else
4941	{
4942	E1=E2;
4943	}
4944	}
4945	kill E2;
4946	list E2;
4947	}
4948	/*
4949	if((!ok)&&(!re)&&(size(#)==1))
4950	{
4951
4952	"the center is wrong";
4953	"it could be one of the following list";
4954	ll;
4955	~;
4956	}
4957	*/
4958	if((!ok)&&(!re)&&(size(#)==2))
4959	{
4960	return(2); //for Center correction
4961	}
4962	return(re);
4963	}
4964	//////////////////////////////////////////////////////////////////////////////
4965	static proc normalCrossB(ideal J,list E,ideal V)
4966	"Internal procedure - no help and no example available
4967	"
4968	{
4969	int i,d,j;
4970	int n=nvars(basering);
4971	list E1,E2;
4972	ideal K,M,Estd;
4973	intvec v,w;
4974
4975	for(i=1;i<=size(E);i++)
4976	{
4977	Estd=std(E[i]+J);
4978	if(deg(Estd[1])>0)
4979	{
4980	E1[size(E1)+1]=Estd;
4981	}
4982	}
4983	E=E1;
4984	for(i=1;i<=size(E);i++)
4985	{
4986	v=i;
4987	E1[i]=list(E[i],v);
4988	}
4989	list ll;
4990	int re=1;
4991
4992	while((size(E1)>0)&&(re==1))
4993	{
4994	K=E1[1][1];
4995	v=E1[1][2];
4996	attrib(K,"isSB",1);
4997	E1=delete(E1,1);
4998	d=n-dim(K);
4999	M=minor(jacob(K),d)+K;
5000	if(deg(std(M+V)[1])>0)
5001	{
5002	re=0;
5003	break;
5004	}
5005	for(i=1;i<=size(E);i++)
5006	{
5007	for(j=1;j<=size(v);j++){if(v[j]==i){break;}}
5008	if(j<=size(v)){if(v[j]==i){i++;continue;}}
5009	Estd=std(K+E[i]);
5010	w=v;
5011	if(deg(Estd[1])==0){i++;continue;}
5012	if(d==n-dim(Estd))
5013	{
5014	if(deg(std(Estd+V)[1])>0)
5015	{
5016	re=0;
5017	break;
5018	}
5019	}
5020	w[size(w)+1]=i;
5021	E2[size(E2)+1]=list(Estd,w);
5022	}
5023	if(size(E2)>0)
5024	{
5025	if(size(E1)>0)
5026	{
5027	E1[size(E1)+1..size(E1)+size(E2)]=E2[1..size(E2)];
5028	}
5029	else
5030	{
5031	E1=E2;
5032	}
5033	}
5034	kill E2;
5035	list E2;
5036	}
5037	return(re);
5038	}
5039	//////////////////////////////////////////////////////////////////////////////
5040	static proc norC(list BO,ideal cent)
5041	"Internal procedure - no help and no example available
5042	"
5043	{
5044	int j;
5045	list E=BO[4];
5046	ideal N=BO[2];
5047	if(BO[3][1]>1){return(0);}
5048	if(deg(std(slocusE(BO[2]))[1])>0){return(0);}
5049	if(!transversalT(N,E)){return(0);}
5050	for(j=1;j<=size(E);j++){if(deg(E[j][1])>0){E[j]=E[j]+N;}}
5051	if(!normalCross(E,cent)){return(0);}
5052	return(1);
5053	}
5054	//////////////////////////////////////////////////////////////////////////////
5055	static proc specialReduce(ideal I,ideal J,poly p)
5056	{
5057	matrix M;
5058	int i,j;
5059	for(i=1;i<=ncols(I);i++)
5060	{
5061	M=coeffs(I[i],@z);
5062	I[i]=0;
5063	for(j=1;j<=nrows(M);j++)
5064	{
5065	I[i]=I[i]+M[j,1]*p^(nrows(M)-j);
5066	}
5067	I[i]=reduce(I[i],J);
5068	}
5069	return(I);
5070	}

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

Download in other formats:

Original Format