Context Navigation

resolve.lib @ fa8122

spielwiese

Last change on this file since fa8122 was fa8122, checked in by Hans Schönemann <hannes@…>, 18 years ago
*hannes: Michael C. git-svn-id: file:///usr/local/Singular/svn/trunk@9349 2c84dea3-7e68-4137-9b89-c4e89433aadc
Property mode set to `100644`
File size: 124.7 KB

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

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

Download in other formats:

Original Format