Context Navigation

derham.lib @ 4c20ee

spielwiese

Last change on this file since 4c20ee was 4c20ee, checked in by Hans Schoenemann <hannes@…>, 11 years ago
fix: new version numbers for libs
Property mode set to `100644`
File size: 121.2 KB

Line
1	////////////////////////////////////////////////////////////////////////////////////
2	version="version derham.lib 4.0.0.0 Jun_2013 ";
3	category="Noncommutative";
4	info="
5	LIBRARY: derham.lib Computation of deRham cohomology
6
7	AUTHORS: Cornelia Rottner, rottner@mathematik.uni-kl.de
8
9	OVERVIEW:
10	A library for computing the de Rham cohomology of complements of complex affine
11	varieties.
12
13
14	REFERENCES:
15	[OT] Oaku, T.; Takayama, N.: Algorithms of D-modules - restriction, tensor product,
16	localzation, and local cohomology groups, J. Pure Appl. Algebra 156, 267-308
17	(2001)
18	[R] Rottner, C.: Computing de Rham Cohomology,diploma thesis (2012)
19	[W1] Walther, U.: Algorithmic computation of local cohomology modules and the local
20	cohomological dimension of algebraic varieties, J. Pure Appl. Algebra 139,
21	303-321 (1999)
22	[W2] Walther, U.: Algorithmic computation of de Rham Cohomology of Complements of
23	Complex Affine Varieties, J. Symbolic Computation 29, 796-839 (2000)
24
25
26	PROCEDURES:
27
28	deRhamCohomology(list[,opt]); computes the de Rham cohomology
29	MVComplex(list); computes the Mayer-Vietoris complex
30	";
31
32	LIB "nctools.lib";
33	LIB "matrix.lib";
34	LIB "qhmoduli.lib";
35	LIB "general.lib";
36	LIB "dmod.lib";
37	LIB "bfun.lib";
38	LIB "dmodapp.lib";
39	LIB "poly.lib";
40	LIB "schreyer.lib";
41
42	////////////////////////////////////////////////////////////////////////////////////
43
44	proc deRhamCohomology(list L,list #)
45	"USAGE: deRhamCohomology(L[,choices]); L a list consisting of polynomials, choices
46	optional list consisting of one up to three strings @*
47	The optional strings may be one of the strings@*
48	-'Vdres': compute the Cartan-Eilenberg resolutions via V_d-homogenization
49	and without using Schreyer's method @*
50	-'Sres': compute the Cartan-Eilenberg resolutions in the homogenized Weyl
51	algebra using Schreyer's method@*
52	one of the strings@*
53	-'iterativeloc': compute localizations by factorizing the polynomials and
54	sucessive localization of the factors @*
55	-'no iterativeloc': compute localizations by directly localizing the
56	product@*
57	and one of the strings
58	-'onlybounds': computes bounds for the minimal and maximal interger roots
59	of the global b-function
60	-'exactroots' computes the minimal and maximal integer root of the global
61	b-function
62	The default is 'Sres', 'iterativeloc' and 'onlybounds'.
63	ASSUME: -The basering must be a polynomial ring over the field of rational numbers@*
64	RETURN: list, where the ith entry is the (i-1)st de Rham cohomology group of the
65	complement of the complex affine variety given by the polynomials in L
66	EXAMPLE:example deRhamCohomology; shows an example
67	"
68	{
69	intvec saveoptions=option(get);
70	intvec i1,i2;
71	option(none);
72	int recursiveloc=1;
73	int i,j,nr,nc;
74	def R=basering;
75	poly islcm, forlcm;
76	int n=nvars(R);
77	int le=size(L)+n;
78	string Syzstring="Sres";
79	int onlybounds=1;
80	for (i=1; i<=size(#); i++)
81	{
82	if (#[i]=="Vdres")
83	{
84	Syzstring="Vdres";
85	}
86	if (#[i]=="noiterativeloc")
87	{
88	recursiveloc=0;
89	}
90	if (#[i]=="exactroots")
91	{
92	onlybounds=0;
93	}
94	}
95	for (i=1; i<=size(L); i++)
96	{
97	if (L[i]==0)
98	{
99	L=delete(L,i);
100	i=i-1;
101	}
102	}
103	if (size(L)==0)
104	{
105	return (list(0));
106	}
107	for (i=1; i<= size(L); i++)
108	{
109	if (leadcoef(L[i])-L[i]==0)
110	{
111	return(list(1));
112	}
113	}
114	if (size(L)==0)
115	{
116	/the complement of the variety given by the input is the whole space/
117	return(list(1));
118	}
119	for (i=1; i<=size(L); i++)
120	{
121	if (typeof(L[i])!="poly")
122	{
123	print("The input list must consist of polynomials");
124	retrun();
125	}
126	}
127	/* 1st step: compute the Mayer-Vietoris Complex and its Fourier transform*/
128	def W=MVComplex(L,recursiveloc);//new ring that contains the MV complex
129	setring W;
130	list fortoVdstrict=MV;
131	ideal IFourier=var(n+1);
132	for (i=2;i<=n;i++)
133	{
134	IFourier=IFourier,var(n+i);
135	}
136	for (i=1; i<=n;i++)
137	{
138	IFourier=IFourier,-var(i);
139	}
140	map cFourier=W,IFourier;
141	matrix sup;
142	for (i=1; i<=size(MV); i++)
143	{
144	sup=fortoVdstrict[i];
145	/takes the Fourier transform of the MV complex/
146	fortoVdstrict[i]=cFourier(sup);
147	}
148	/* 2nd step: Compute a V_d-strict free complex that is quasi-isomorphic to the
149	complex fortoVdstrict
150	The 1st entry of the list rem will be the quasi-isomorphic complex, the 2nd
151	entry contains the cohomology modules and is needed for the computation of the
152	global b-function*/
153	list rem=toVdStrictFreeComplex(fortoVdstrict,Syzstring);
154	list newcomplex=rem[1];
155	////////////////////////////////////////////////////////////////////////////////////
156	/* 3rd step: Compute the bounds for the minimal and maximal integer root of the
157	global b-function of newcomplex(i.e. compute the lcm of the b-functions of its
158	cohomology modules)(if onlybouns=1). Else we compute the minimal and maximal
159	integer root.
160
161	If we compute only the bounds, we omit additional Groebner basis computations.
162	However this leads to a higher-dimensional truncated complex.
163
164	Note that the cohomology modules are already contained in rem[2].
165	minmaxk[1] and minmaxk[2] will contain the bounds resp exact roots.*/
166	if (onlybounds==1)
167	{
168	list minmaxk=globalBFun(rem[2],Syzstring);
169	}
170	else
171	{
172	list minmaxk=exactGlobalBFun(rem[2],Syzstring);
173	}
174	if (size(minmaxk)==0)
175	{
176	return (0);
177	}
178	/*4th step: Truncate the complex D_n/(x_1,...,x_n)\otimes C, (where
179	C=(C^i[m^i],d^i) is given by newcomplex, i.e. C^i=D_n^newcomplex[3*i-2],
180	m^i=newcomplex[3i-1], d^i=newcomplex[3i]), using Thm 5.7 in [W1]:
181	The truncated module D_n/(x_1,..,x_n)\otimes C[i] is generated by the set
182	(0,...,P_(i_j),0,...), where P_(i_j) is a monomial in C[D(1),...,D(n)] and
183	if it is placed in component k it holds that
184	minmaxk[1]-m^i[k]<=deg(P_(i_j))<=minmaxk[2]-m^i[k]*/
185	int k,l;
186	list truncatedcomplex,shorten,upto;
187	for (i=1; i<=size(newcomplex) div 3; i++)
188	{
189	shorten[3*i-1]=list();
190	for (j=1; j<=size(newcomplex[3*i-1]); j++)
191	{
192	/shorten[3i-1][j][k]=minmaxk[k]-m^i[j]+1 (for k=1,2) if this value is
193	positive otherwise we will set it to be list();
194	we added +1, because we will use a list, where we put in position l
195	polys of degree l+1*/
196	shorten[3i-1][j]=list(minmaxk[1]-newcomplex[3i-1][j]+1);
197	shorten[3i-1][j][2]=minmaxk[2]-newcomplex[3i-1][j]+1;
198	upto[size(upto)+1]=shorten[3*i-1][j][2];
199	if (shorten[3*i-1][j][2]<=0)
200	{
201	shorten[3*i-1][j]=list();
202	}
203	else
204	{
205	if (shorten[3*i-1][j][1]<=0)
206	{
207	shorten[3*i-1][j][1]=1;
208	}
209	}
210	}
211	}
212	int iupto=Max(upto);//maximal degree +1 of the polynomials we have to consider
213	if (iupto<=0)
214	{
215	return(list(0));
216	}
217	list allpolys;
218	/allpolys[i] will consist list of all monomials in D(1),...,D(n) of degree i-1/
219	allpolys[1]=list(1);
220	list minvar;
221	minvar[1]=list(1);
222	for (i=1; i<=iupto-1; i++)
223	{
224	allpolys[i+1]=list();
225	minvar[i+1]=list();
226	for (k=1; k<=size(allpolys[i]); k++)
227	{
228	for (j=minvar[i][k]; j<=nvars(W) div 2; j++)
229	{
230	allpolys[i+1][size(allpolys[i+1])+1]=allpolys[i][k]*D(j);
231	minvar[i+1][size(minvar[i+1])+1]=j;
232	}
233	}
234	}
235	list keepformatrix,sizetruncom,fortrun,fst;
236	int count,stc;
237	intvec v,forin;
238	matrix subm;
239	/now we compute the truncation/
240	for (i=1; i<=size(newcomplex) div 3; i++)
241	{
242	/truncatedcomplex[2i-1] will contain all the generators for the truncation
243	of D_n/(x(1),..,x(n))\otimes C[i]*/
244	truncatedcomplex[2*i-1]=list();
245	sizetruncom[2*i-1]=list();
246	sizetruncom[2*i]=list();
247	/truncatedcomplex[2i] will be the map trunc(D_n/(x(1),..,x(n))\otimes C[i])
248	->trunc(D_n/(x(1),..,x(n))\otimes C[i+1])*/
249	truncatedcomplex[2i]=newcomplex[3i];
250	v=0;count=0;
251	sizetruncom[2*i][1]=0;
252	for (j=1; j<=newcomplex[3*i-2]; j++)
253	{
254	if (size(shorten[3*i-1][j])!=0)
255	{
256	fortrun=sublist(allpolys,shorten[3i-1][j][1],shorten[3i-1][j][2]);
257	truncatedcomplex[2i-1][size(truncatedcomplex[2i-1])+1]=fortrun[1];
258	count=count+fortrun[2];
259	fst=list(int(shorten[3i-1][j][1])-1,int(shorten[3i-1][j][2])-1);
260	sizetruncom[2i-1][size(sizetruncom[2i-1])+1]=fst;
261	sizetruncom[2i][size(sizetruncom[2i])+1]=count;
262	if (v!=0)
263	{
264	v[size(v)+1]=j;
265	}
266	else
267	{
268	v[1]=j;
269	}
270	}
271	}
272	if (v!=0)
273	{
274	subm=submat(truncatedcomplex[2i],v,1..ncols(truncatedcomplex[2i]));
275	truncatedcomplex[2*i]=subm;
276	if (i!=1)
277	{
278	i1=1..nrows(truncatedcomplex[2*(i-1)]);
279	subm=submat(truncatedcomplex[2*(i-1)],i1,v);
280	truncatedcomplex[2*(i-1)]=subm;
281	}
282	}
283	else
284	{
285	truncatedcomplex[2i]=matrix(0,1,ncols(truncatedcomplex[2i]));
286	if (i!=1)
287	{
288	nr=nrows(truncatedcomplex[2*(i-1)]);
289	truncatedcomplex[2*(i-1)]=matrix(0,nr,1);
290	}
291	}
292	}
293	matrix M;
294	int st,pi,pj;
295	poly ptc;
296	int b,d,ideg,kplus,lplus;
297	poly form,lform,nform;
298	/computation of the maps/
299	for (i=1; i<size(truncatedcomplex) div 2; i++)
300	{
301	nr=max(1,sizetruncom[2i][size(sizetruncom[2i])]);
302	nc=max(1,sizetruncom[2i+2][size(sizetruncom[2i+2])]);
303	M=matrix(0,nr,nc);
304	for (k=1; k<=size(truncatedcomplex[2*i-1]);k++)
305	{
306	for (l=1; l<=size(truncatedcomplex[2*(i+1)-1]); l++)
307	{
308	if (size(sizetruncom[2*i])!=1)
309	{
310	for (j=1; j<=size(truncatedcomplex[2*i-1][k]); j++)
311	{
312	for (b=1; b<=size(truncatedcomplex[2*i-1][k][j]); b++)
313	{
314	form=truncatedcomplex[2*i-1][k][j][b][1];
315	form=formtruncatedcomplex[2i][k,l];
316	while (form!=0)
317	{
318	lform=lead(form);
319	v=leadexp(lform);
320	v=v[1..n];
321	if (v==(0:n))
322	{
323	ideg=deg(lform)-sizetruncom[2*(i+1)-1][l][1];
324	if (ideg>=0)
325	{
326	nr=ideg+1;
327	st=size(truncatedcomplex[2*(i+1)-1][l][nr]);
328	for (d=1; d<=st;d++)
329	{
330	nc=2*(i+1)-1;
331	ptc=truncatedcomplex[nc][l][ideg+1][d][1];
332	if (leadmonom(lform)==ptc)
333	{
334	nr=2*i-1;
335	pi=truncatedcomplex[nr][k][j][b][2];
336	pi=pi+sizetruncom[2*i][k];
337	nc=2*(i+1)-1;
338	nr=ideg+1;
339	pj=truncatedcomplex[nc][l][nr][d][2];
340	pj=pj+sizetruncom[2*(i+1)][l];
341	M[pi,pj]=leadcoef(lform);
342	break;
343	}
344	}
345	}
346	}
347	form=form-lform;
348	}
349	}
350	}
351	}
352	}
353	}
354	truncatedcomplex[2*i]=M;
355	truncatedcomplex[2i-1]=sizetruncom[2i][size(sizetruncom[2*i])];
356	}
357	truncatedcomplex[2i-1]=sizetruncom[2i][size(sizetruncom[2*i])];
358	if (truncatedcomplex[2*i-1]!=0)
359	{
360	truncatedcomplex[2i]=matrix(0,truncatedcomplex[2i-1],1);
361	}
362	setring R;
363	list truncatedcomplex=imap(W,truncatedcomplex);
364	/*computes the cohomology of the complex (D^i,d^i) given by truncatedcomplex,
365	i.e. D^i=C^truncatedcomplex[2i-1] and d^i=truncatedcomplex[2i]*/
366	list derhamhom=findCohomology(truncatedcomplex,le);
367	option(set,saveoptions);
368	return (derhamhom);
369	}
370
371	example
372	{ "EXAMPLE:";
373	ring r = 0,(x,y,z),dp;
374	list L=(xy,xz);
375	deRhamCohomology(L);
376	}
377
378	////////////////////////////////////////////////////////////////////////////////////
379	//COMPUTATION OF THE MAYER-VIETORIS COMPLEX
380	////////////////////////////////////////////////////////////////////////////////////
381
382	proc MVComplex(list L,list #)
383	"USAGE:MVComplex(L); L a list of polynomials
384	ASSUME: -Basering is a polynomial ring with n vwariables and rational coefficients
385	-L is a list of non-constant polynomials
386	RETURN: ring W: the nth Weyl algebra @*
387	W contains a list MV, which represents the Mayer-Vietrois complex (C^i,d^i) of the
388	polynomials contained in L as follows:@*
389	the C^i are given by D_n^ncols(C[2i-1])/im(C[2i-1]) and the differentials
390	d^i are given by C[2*i]
391	EXAMPLE:example MVComplex; shows an example
392	"
393	{
394	/* We follow algorithm 3.2.5 in [R],if #!=0 we use also Remark 3.2.6 in [R] for
395	an additional iterative localization*/
396	def R=basering;
397	int i;
398	int iterative=1;
399	if (size(#)!=0)
400	{
401	iterative=#[1];
402	}
403	for (i=1; i<=size(L); i++)
404	{
405	if (L[i]==0)
406	{
407	print("localization with respect to 0 not possible");
408	return();
409	}
410	if (leadcoef(L[i])-L[i]==0)
411	{
412	print("polynomials must be non-constant");
413	return();
414	}
415	}
416	if (iterative==1)
417	{
418	/*compute the localizations by factorizing the polynomials and iterative
419	localization of the factors */
420	for (i=1; i<=size(L); i++)
421	{
422	L[i]=factorize(L[i],1);
423	}
424	}
425	int r=size(L);
426	int n=nvars(basering);
427	int le=size(L)+n;
428	/construct the ring Ws/
429	def W=makeWeyl(n);
430	setring W;
431	list man=ringlist(W);
432	if (n==1)
433	{
434	man[2][1]="x(1)";
435	man[2][2]="D(1)";
436	def Wi=ring(man);
437	setring Wi;
438	kill W;
439	def W=Wi;
440	setring W;
441	list man=ringlist(W);
442	}
443	man[2][size(man[2])+1]="s";;
444	man[3][3]=man[3][2];
445	man[3][2]=list("dp",intvec(1));
446	matrix N=UpOneMatrix(size(man[2]));
447	man[5]=N;
448	matrix M[1][1];
449	man[6]=transpose(concat(transpose(concat(man[6],M)),M));
450	def Ws=ring(man);
451	setring Ws;
452	int j,k,l,c;
453	list L=fetch(R,L);
454	list Cech;
455	ideal J=var(1+n);
456	for (i=2; i<=n; i++)
457	{
458	J=J,var(i+n);
459	}
460	Cech[1]=list(J);
461	list Theta, remminroots;
462	Theta[1]=list(list(list(),1,1));
463	list rem,findminintroot,diffmaps;
464	int minroot,st,sk;
465	intvec k1;
466	poly fred,forfetch;
467	matrix subm;
468	int rmr;
469	if (iterative==0)
470	{/computation of the modules of the MV complex/
471	for (i=1; i<=r; i++)
472	{
473	findminintroot=list();
474	Cech[i+1]=list();
475	Theta[i+1]=list();
476	k1=1;
477	for (j=1; j<=i; j++)
478	{
479	k1[size(k1)+1]=size(Theta[j+1]);
480	for (k=1; k<=k1[j]; k++)
481	{
482	Theta[j+1][size(Theta[j+1])+1]=list(Theta[j][k][1]+list(i));
483	Theta[j+1][size(Theta[j+1])][2]=Theta[j][k][2]*L[i];
484	/*We compute the s-parametric annihilator J(s) and the b-function
485	of the polynomial L[i] and Cech[i][k] to localize the module
486	D_n/(D(1),...,D(n))[L[i]^(-1)]\otimes D_n^c/im(Cech[i][k]),
487	where c=ncols(Cech[i][k]) and the im(Cech[i][k]) is generated by
488	the rows of the matrix.
489	If we plug the minimal integer root r(or a smaller integer
490	value)in J(s), then D_n^ncols(J(s))/im(J(r)) is isomorphic to
491	the above localization*/
492	rem=SannfsIBM(L[i],Cech[j][k]);
493	Cech[j+1][size(Cech[j+1])+1]=rem[1];
494	findminintroot[size(findminintroot)+1]=rem[2];
495	}
496	}
497	/* we compute the minimal root of all b-functions of L[i] computed above,
498	because we want to plug in the same root r in all s-parametric
499	annihilators we computed for L[i] ->this will ensure we can compute
500	the maps of the MV complex*/
501	minroot=minIntRoot0(findminintroot);
502	for (j=1; j<=i; j++)
503	{
504	for (k=1; k<=k1[j]; k++)
505	{
506	sk=size(Cech[j+1])+1-k;
507	Cech[j+1][size(Cech[j+1])+1-k]=subst(Cech[j+1][sk],s,minroot);
508	}
509	}
510	remminroots[i]=minroot;
511	}
512	Cech=delete(Cech,1);
513	Theta=delete(Theta,1);
514	list zw;
515	poly reme;
516	/computation of the maps of the MV complex/
517	for (i=1; i<r; i++)
518	{
519	diffmaps[i]=matrix(0,size(Cech[i]),size(Cech[i+1]));
520	for (j=1; j<=size(Cech[i]); j++)
521	{
522	for (k=1; k<=size(Cech[i+1]); k++)
523	{
524	zw=LMSubset(Theta[i][j][1],Theta[i+1][k][1]);
525	if (zw[2]!=0)
526	{
527	rmr=-remminroots[zw[1]];
528	reme=zw[2]*(Theta[i+1][k][2]/Theta[i][j][2])^(rmr);
529	zw[2]=zw[2]*(Theta[i+1][k][2]/Theta[i][j][2])^(rmr);
530	diffmaps[i][j,k]=zw[2];
531	}
532	}
533	}
534	}
535	diffmaps[r]=matrix(0,1,1);
536	}
537	if (iterative==1)
538	{
539	for (i=1; i<=r;i++)
540	{
541	Cech[i+1]=list();
542	Theta[i+1]=list();
543	k1=1;
544	for (c=1; c<=size(L[i]); c++)
545	{
546	findminintroot=list();
547	for (j=1; j<=i; j++)
548	{
549	if (c==1)
550	{
551	k1[size(k1)+1]=size(Theta[j+1]);
552	}
553	for (k=1; k<=k1[j]; k++)
554	{
555	/*We compute the s-parametric annihilator J(s) und the b-
556	function of the polynomial L[i][c] and Cech[i][k] to
557	localize the module D_n/(D(1),...,D(n))[L[i][c]^(-1)]\otimes
558	D_n^c/im(Cech[i][k]), where c=ncols(Cech[i][k]).
559	If we plug the minimal integer root r(or a smaller integer
560	value)in J(s), then D_n^ncols(J(s))/im(J(r)) is isomorphic
561	to the above localization*/
562	if (c==1)
563	{
564	rmr=size(Theta[j+1])+1;
565	Theta[j+1][rmr]=list(Theta[j][k][1]+list(i));
566	Theta[j+1][size(Theta[j+1])][2]=Theta[j][k][2]*L[i][c];
567	rem=SannfsIBM(L[i][c],Cech[j][k]);
568	Cech[j+1][size(Cech[j+1])+1]=rem[1];
569	findminintroot[size(findminintroot)+1]=rem[2];
570	}
571	else
572	{
573	st=size(Theta[j+1])-k1[j]+k;
574	Theta[j+1][st][2]=Theta[j+1][st][2]*L[i][c];
575	rem=SannfsIBM(L[i][c],Cech[j+1][size(Cech[j+1])-k1[j]+k]);
576	Cech[j+1][size(Cech[j+1])-k1[j]+k]=rem[1];
577	findminintroot[size(findminintroot)+1]=rem[2];
578	}
579	}
580	}
581	/* we compute the minimal root of all b-functions of L[i][c]
582	computed above,because we want to plug in the same root r in all
583	s-parametric annihilators we computed for L[i] ->this will
584	ensure we can compute the maps of the MV complex*/
585	minroot=minIntRoot0(findminintroot);
586	for (j=1; j<=i; j++)
587	{
588	for (k=1; k<=k1[j]; k++)
589	{
590	st=size(Cech[j+1])+1-k;
591	Cech[j+1][st]=subst(Cech[j+1][st],s,minroot);
592	}
593	}
594	if (c==1)
595	{
596	remminroots[i]=list();
597	}
598	remminroots[i][c]=minroot;
599	}
600	}
601	Cech=delete(Cech,1);
602	Theta=delete(Theta,1);
603	list zw;
604	poly reme;
605	/maps of the MV Complex/
606	for (i=1; i<r; i++)
607	{
608	diffmaps[i]=matrix(0,size(Cech[i]),size(Cech[i+1]));
609	for (j=1; j<=size(Cech[i]); j++)
610	{
611	for (k=1; k<=size(Cech[i+1]); k++)
612	{
613	zw=LMSubset(Theta[i][j][1],Theta[i+1][k][1]);
614	if (zw[2]!=0)
615	{
616	reme=1;
617	for (c=1; c<=size(L[zw[1]]);c++)
618	{
619	reme=reme*L[zw[1]][c]^(-remminroots[zw[1]][c]);
620	}
621	diffmaps[i][j,k]=zw[2]*reme;
622	}
623	}
624	}
625	}
626	diffmaps[r]=matrix(0,1,1);
627	}
628	setring W;
629	/map the modules and maps to the Weyl algebra/
630	list diffmaps=imap(Ws,diffmaps);
631	list Cechmodules=imap(Ws,Cech);
632	list Cech;
633	matrix sup;
634	for (i=1; i<=r; i++)
635	{
636	sup=transpose(matrix(Cechmodules[i][1]));
637	Cech[2*i-1]=sup;
638	for (j=2; j<=size(Cechmodules[i]); j++)
639	{
640	sup=transpose(matrix(Cechmodules[i][j]));
641	Cech[2i-1]=dsum(Cech[2i-1],sup);
642	}
643	sup=matrix(diffmaps[i]);
644	Cech[2*i]=sup;
645	}
646	list MV=Cech;
647	export MV;
648	return (W);
649	}
650
651	example
652	{ "EXAMPLE:";
653	ring r = 0,(x,y,z),dp;
654	list L=xy,xz;
655	def C=MVComplex(L);
656	setring C;
657	MV;
658	}
659
660	////////////////////////////////////////////////////////////////////////////////////
661
662	static proc SannfsIBM(poly F,ideal myJ)
663	"USAGE: SannfsIBM(f,J), F poly, J ideal
664	ASSUME: basering is D_n[s], where D_n is the Weyl algebra and s and extra
665	commutative variable@*
666	f is a polynomial in the variables x(1),...,x(n) with rational coefficients
667	@*
668	J is holonomic and f-saturated
669	RETURN AlList of the form (K,g), where K is an ideal and g a univariant polynomial
670	in the variable s. K is the s-parametric annihilator of F and J and g is
671	the b-function of F and J.
672	"
673	{
674	/*modified version of the procedure SannfsBM from the library dmod.lib: SannfsBM
675	computes the s-parametric annihilator for J=(x_1,...,x_n)*/
676	/* We use Algorithm 3.1.12 in[R] to compute the s-parametric
677	annihilator. Then we use the s-parametric annihilator to compute the b-function
678	via Algorithm 4.7 in [W1].*/
679	/* We assume that the basering the the nth Weyl algebra D_n. We create the ring
680	D_n[s,t], where ts=st-t*/
681	def save = basering;
682	int N = nvars(basering)-1;
683	int Nnew = N+2;
684	int i,j;
685	string s;
686	list RL = ringlist(basering);
687	list L, Lord;
688	list tmp;
689	intvec iv;
690	L[1] = RL[1];
691	L[4] = RL[4];
692	list Name = RL[2];
693	Name=delete(Name,size(Name));
694	list RName;
695	RName[1] = "t";
696	RName[2] = "s";
697	list DName;
698	for(i=1;i<=N div 2;i++)
699	{
700	DName[i] = var(N div 2+i);
701	Name=delete(Name,N div 2+1);
702	}
703	tmp[1] = "t";
704	tmp[2] = "s";
705	list NName = tmp +Name+DName;
706	L[2] = NName;
707	kill NName;
708	tmp[1] = "lp";
709	iv = 1,1;
710	tmp[2] = iv;
711	Lord[1] = tmp;
712	tmp[1] = "dp";
713	s = "iv=";
714	for(i=1;i<=Nnew;i++)
715	{
716	s = s+"1,";
717	}
718	s[size(s)]= ";";
719	execute(s);
720	kill s;
721	tmp[2] = iv;
722	Lord[2] = tmp;
723	tmp[1] = "C";
724	iv = 0;
725	tmp[2] = iv;
726	Lord[3] = tmp;
727	tmp = 0;
728	L[3] = Lord;
729	def @R@ = ring(L);
730	setring @R@;
731	matrix @D[Nnew][Nnew];
732	@D[1,2]=t;
733	for(i=1; i<=N div 2; i++)
734	{
735	@D[2+i, N div 2+2+i]=1;
736	}
737	def @R = nc_algebra(1,@D);
738	setring @R;
739	kill @R@;
740	/we start with the computation of the s-parametric annihilator/
741	poly F = imap(save,F);
742	ideal myJ=imap(save,myJ);
743	for (i=1; i<=N div 2; i++)
744	{
745	myJ=subst(myJ,D(i),D(i)+diff(F,x(i))*t);
746	}
747	ideal I = t*F+s;
748	I=I,myJ;//the s-parametric annihilator in D_n[s,t]
749	/we compute the intersection of I and D_n[s]/
750	ideal J = slimgb(I);
751	ideal K = nselect(J,1);
752	K = slimgb(K);//the s-parametric annihilator
753	/we use K to compute the b-function/
754	ideal B=K,F;
755	B=slimgb(B);
756	vector p=pIntersect(s,B);
757	poly f=vec2poly(p,2);
758	setring save;
759	poly f=imap(@R,f);
760	ideal K=imap(@R,K);
761	return (list(K,f));
762	}
763
764	////////////////////////////////////////////////////////////////////////////////////
765	//COMPUTATION OF QA UASI-ISOMORPHIC V_D-STRICT FREE COMPLEX
766	////////////////////////////////////////////////////////////////////////////////////
767
768	static proc toVdStrictFreeComplex(list L,string Syzstring,list #)
769	"USAGE: toVdStrictFreeComplex(L, Syzstring [,d]); L a list of the form
770	(M_1,f_1,...,M_s,f_s), where the M_i and f_i are matrices, Syzstring a
771	string, d an optional integer
772	ASSUME: Basering is the Weyl algebra D_n @*
773	(M_1,f_1,...,M_s,f_s) represents a complex 0->D_n^(r_1)/im(M_1)->
774	D_n^(r_2)/im(M_2)->...->D_n^(r_s)->0 with differentials f_i, where im(M_i)
775	is generated by the rows of M_i. In particular it hold:@*
776	- The M_i are m_i x r_i-matrices and the f_iare r_i x r_(i+1)-matrices @*
777	-the image of M_1f_i is contained in the image of M_(i+1) @
778	d is an integer between 1 and n. If no value for d is given, it is assumed
779	to be n @*
780	Syzstring is either: @*
781	-'Sres' (computes the resolutions and Groebner bases in the homogenized
782	Weyl algebra using Schreyer's method)@*
783	or @*
784	-'Vdres' (computes the resolutions via V_d-homogenization and without
785	Schreyer's method)@*
786	RETURN: list of the form (L_1,L_2), were L_1 and L_2 are lists @*
787	L_1 is of the form (i_(-n-1),g_(-n-1),m_(-n-1),...,i_s,g_s,m_s) such that:@*
788	-the i_j are integers, the g_j are i_j x i_(j+1)-matrices, the m_j intvecs
789	of size i_j@*
790	-D_n^(i_(-n-1))[m_(-n-1)]->...->D_n^(i_s)[m_s]->0 is a V_d-strict complex
791	with differentials m_i that is quasi-isomorphic to the complex given by L@*
792	L_2 is of the form (H_1,n_1,...,H_s,n_s), where the H_i are matrices and
793	the n_i are shift vectors such that:@*
794	-coker(H_i) is the ith cohomology group of the complex given by L_1@*
795	-the n_i are the shift vectors of the coker(H_i)
796	THEORY: We follow Algorithm 3.8 in [W2]
797	"
798	{
799	def B=basering;
800	int n=nvars(B) div 2+2;
801	int d=nvars(B) div 2;
802	intvec v;
803	list out, outall;
804	int i,j,k,indi,nc,nr;
805	matrix mem;
806	intvec i1,i2;
807	if (size(#)!=0)
808	{
809	for (i=1; i<=size(#); i++)
810	{
811	if (typeof(#[i])=="int")
812	{
813	if (#[1]>=1 and #[1]<=n)
814	{
815	d=#[i];
816	}
817	}
818	}
819	}
820	/* If size(L)=2, our complex consists for only one non-trivial module.
821	Therefore, we just have to compute a V_d-strict resolution of this module.*/
822	if (size(L)==2)
823	{
824	v=(0:ncols(L[1]));
825	out[3*n-1]=v;
826	out[3*n-2]=ncols(L[1]);
827	out[3*n]=L[2];
828	if (Syzstring=="Vdres")
829	{
830	/*if Syzstring="Vdres", we compute a V_d-strict Groebner basis of L[1]
831	using F-homogenization (Prop. 3.9 in [OT]); then we compute the syzygies
832	and make them V_d-strict using Prop 3.9[OT] and so on*/
833	out[3*n-3]=VdStrictGB(L[1],d,v);
834	for (i=n-1; i>=1; i--)
835	{
836	out[3i-2]=nrows(out[3i]);
837	v=0;
838	for (j=1; j<=out[3*i-2]; j++)
839	{
840	mem=submat(out[3i],j,intvec(1..ncols(out[3i])));
841	v[j]=VdDeg(mem,d, out[3*i+2]);//next shift vector
842	}
843	out[3*i-1]=v;
844	if (i!=1)
845	{
846	/next step in the resolution/
847	out[3i-3]=transpose(syz(transpose(out[3i])));
848	if (out[3i-3]!=matrix(0,nrows(out[3i-3]),ncols(out[3*i-3])))
849	{
850	/makes the resolution V_d-strict/
851	out[3i-3]=VdStrictGB(out[3i-3],d,out[3*i-1]);
852	}
853	else
854	{
855	/resolution is already computed/
856	out[3i-3]=matrix(0,1,ncols(out[3i-3]));
857	out[3*i-4]=intvec(0);
858	out[3*i-5]=int(0);
859	for (j=i-2; j>=1; j--)
860	{
861	out[3*j]=matrix(0,1,1);
862	out[3*j-1]=intvec(0);
863	out[3*j-2]=int(0);
864	}
865	break;
866	}
867	}
868	}
869	}
870	else
871	{
872	/*in the case Syzstring!="Vdres" we compute the resolution in the
873	homogenized Weyl algebra using Thm 9.10 in[OT]*/
874	def HomWeyl=makeHomogenizedWeyl(d);
875	setring HomWeyl;
876	list L=fetch(B,L);
877	L[1]=nHomogenize(L[1]);
878	list out=fetch(B,out);
879	out[3*n-3]=L[1];
880	/*computes a ring with a list RES; RES is a V_d-strict resolution of
881	L[1]*/
882	def ringofSyz=Sres(transpose(L[1]),d);
883	setring ringofSyz;
884	int logens=2;
885	matrix mem;
886	list out=fetch(HomWeyl,out);
887	out[3*n-3]=transpose(matrix(RES[2]));
888	out[3n-3]=subst(out[3n-3],h,1);
889	for (i=n-1; i>=1; i--)
890	{
891	out[3i-2]=nrows(out[3i]);
892	v=0;
893	for (j=1; j<=out[3*i-2]; j++)
894	{
895	mem=submat(out[3i],j,intvec(1..ncols(out[3i])));
896	v[j]=VdDeg(mem,d, out[3*i+2]);
897	}
898	out[3*i-1]=v;//shift vector such that the resolution RES is V_d-strict
899	if (i!=1)
900	{
901	indi=0;
902	if (size(RES)>=n-i+2)
903	{
904	nr=nrows(matrix(RES[n-i+2]));
905	mem=matrix(0,nr,ncols(matrix(RES[n-i+2])));
906	if (matrix(RES[n-i+2])!=mem)
907	{
908	indi=1;
909	out[3*i-3]=(matrix(RES[n-i+2]));
910	if (nrows(out[3i-3])-logens+1!=nrows(out[3i]))
911	{
912	mem=out[3*i-3];
913	out[3*i-3]=matrix(mem,nrows(mem)+logens-1,ncols(mem));
914	}
915	mem=out[3*i-3];
916	i1=intvec(logens..nrows(mem));
917	mem=submat(mem,i1,intvec(1..ncols(mem)));
918	out[3*i-3]=transpose(mem);
919	out[3i-3]=subst(out[3i-3],h,1);
920	logens=logens+ncols(out[3*i-3]);
921	}
922	}
923	if(indi==0)
924	{
925	out[3i-3]=matrix(0,1,nrows(out[3i]));
926	out[3*i-4]=intvec(0);
927	out[3*i-5]=int(0);
928	for (j=i-2; j>=1; j--)
929	{
930	out[3*j]=matrix(0,1,1);
931	out[3*j-1]=intvec(0);
932	out[3*j-2]=int(0);
933	}
934	break;
935	}
936	}
937	}
938	setring B;
939	out=fetch(ringofSyz,out);//contains the V_d-strict resolution
940	kill ringofSyz;
941	}
942	outall[1]=out;
943	outall[2]=list(list(out[3n-3],out[3n-1]));
944	return(outall);
945	}
946	/*case size(L)>2: We compute a quasi-isomorphic free complex following Alg 3.8 in
947	[W2]*/
948	/* We denote the complex given by L as (C^i,d^i).
949	We start by computing in the proc shortExaxtPieces representations for the
950	short exact sequences B^i->Z^i->H^i and Z^i->C^i->B^(i+1), where the B^i, Z^i
951	and H^i are coboundaries, cocycles and cohomology groups, respectively.*/
952	out=shortExactPieces(L);
953	list rem;
954	/* shortExactpiecesToVdStrict makes the sequences B^i->Z^i->H^i and
955	Z^i->C^i->B^(i+1) V_d-strict*/
956	rem=shortExactPiecesToVdStrict(out,d,Syzstring);
957	/*VdStrictDoubleComplexes computes V_d-strict resolutions over the seqeunces from
958	proc shortExactPiecesToVdstrict*/
959	out=VdStrictDoubleComplexes(rem[1],d,Syzstring);
960	for (i=1;i<=size(out); i++)
961	{
962	rem[2][i][1]=out[i][1][5][1];
963	rem[2][i][2]=out[i][1][8][1];
964	}
965	/* AssemblingDoubleComplexes puts the resolution of the C^i (from the sequences
966	Z^i->C^i->B^(i+1)) together to obtain a Cartan-Eilenberg resolution of
967	(C^i,d^i)*/
968	out=assemblingDoubleComplexes(out);
969	/*the proc totalComplex takes the total complex of the double complex from the
970	proc assemblingDoubleComplexes*/
971	out=totalComplex(out);
972	outall[1]=out;
973	outall[2]=rem[2];//contains the cohomology groups and their shift vectors
974	return (outall);
975	}
976
977	////////////////////////////////////////////////////////////////////////////////////
978
979
980	static proc sublist(list L,int m,int n)
981	{
982	list out;
983	int i; int j;
984	int count;
985	for (i=m; i<=n; i++)
986	{
987	out[size(out)+1]=list();
988	for (j=1; j<=size(L[i]); j++)
989	{
990	count=count+1;
991	out[size(out)][j]=list(L[i][j],count);
992	}
993	}
994	list o=list(out,count);
995	return(o);
996	}
997
998	////////////////////////////////////////////////////////////////////////////////////
999
1000	static proc LMSubset(list L,list M)
1001	{
1002	int i;
1003	int j=1;
1004	list position=(M[size(M)],(-1)^(size(L)));
1005	for (i=1; i<=size(L); i++)
1006	{
1007	if (L[i]!=M[j])
1008	{
1009	if (L[i]!=M[i+1] or j!=i)
1010	{
1011	return (L[i],0);
1012	}
1013	else
1014	{
1015	position=(M[i],(-1)^(i-1));
1016	j=j+1;
1017	}
1018	}
1019	j=j+1;
1020
1021	}
1022	return (position);
1023	}
1024
1025	////////////////////////////////////////////////////////////////////////////////////
1026
1027	static proc shortExactPieces(list L)
1028	{
1029	/we follow Section 3.3 in [W2]/
1030	/* we assume that L=(M_1,f_1,...,M_s,f_s) defines the complex C=(C^i,d^i)
1031	as in the procedure toVdstrictcomplex*/
1032	matrix Bnew= divdr(L[2],L[3]);
1033	matrix Bold=Bnew;
1034	matrix Z=divdr(Bnew,L[1]);
1035	list bzh,zcb;
1036	bzh=list(list(),list(),Z,unitmat(ncols(Z)),Z);
1037	zcb=(Z, Bnew, L[1], unitmat(ncols(L[1])), Bnew);
1038	list sep;
1039	/* the list sep will be of size s such that
1040	-sep[i]=(sep[i][1],sep[i][2]) is a list of two lists
1041	-sep[i][1]=(B^i,f^(BZi),Z^i,f_^(ZHi),H^i) such that coker(B^i)->coker(Z^i)
1042	->coker(H^i) represents the short exact seqeuence B^i(C)->Z^i(C)->H^i(C)
1043	-sep[i][2]=(Z^i,f^(ZCi),C^i,f^(CBi),B^(i+1)) such that coker(Z^i)->coker(C^i)->
1044	coker(B^(i+1)) represents the short exact seqeuence Z^i(C)->C^i->B^(i+1)(C)*/
1045	sep[1]=list(bzh,zcb);
1046	int i;
1047	list out;
1048	for (i=3; i<=size(L)-2; i=i+2)
1049	{
1050	/the proc bzhzcb computes representations for the short exact seqeunces /
1051	out=bzhzcb(Bold, L[i-1] , L[i], L[i+1], L[i+2]);
1052	sep[size(sep)+1]=out[1];
1053	Bold=out[2];
1054	}
1055	bzh=(divdr(L[size(L)-2], L[size(L)-1]),L[size(L)-2], L[size(L)-1]);
1056	bzh[4]=unitmat(ncols(L[size(L)-1]));
1057	bzh[5]=transpose(concat(transpose(L[size(L)-2]),transpose(L[size(L)-1])));
1058	zcb=(L[size(L)-1], unitmat(ncols(L[size(L)-1])), L[size(L)-1],list(),list());
1059	sep[size(sep)+1]=list(bzh,zcb);
1060	return(sep);
1061	}
1062
1063	////////////////////////////////////////////////////////////////////////////////////
1064
1065	static proc bzhzcb (matrix Bold,matrix f0,matrix C1,matrix f1,matrix C2)
1066	{
1067	matrix Bnew=divdr(f1,C2);
1068	matrix Z= divdr(Bnew,C1);
1069	matrix lift1= matrixLift(Bnew,f0);
1070	matrix H=transpose(concat(transpose(lift1),transpose(Z)));
1071	list bzh=(Bold, lift1, Z, unitmat(ncols(Z)),H);
1072	list zcb=(Z, Bnew, C1, unitmat(ncols(C1)),Bnew);
1073	list out=(list(bzh, zcb), Bnew);
1074	return(out);
1075	}
1076
1077	////////////////////////////////////////////////////////////////////////////////////
1078
1079	static proc shortExactPiecesToVdStrict(list C,int d,list #)
1080	{/* We transform the short exact pieces from procedure shortExactPieces to V_d-
1081	strict short exact sequences. For this, we use Algorithm 3.11 and Lemma 4.2 in
1082	[W2].*/
1083	/* If we compute our Groebner bases in the homogenized Weyl algebra, we already
1084	compute some resolutions it omit additional Groebner basis computations later
1085	on.*/
1086	int s =size(C);int i; int j;
1087	string Syzstring="Sres";
1088	intvec v=0:ncols(C[s][1][5]);
1089	if (size(#)!=0)
1090	{
1091	for (i=1; i<=size(#); i++)
1092	{
1093	if (typeof(#[i])=="string")
1094	{
1095	Syzstring=#[i];
1096	}
1097	if (typeof(#[i])=="intvec")
1098	{
1099	v=#[i];
1100	}
1101	}
1102	}
1103	list out;
1104	list forout;
1105	if (Syzstring=="Vdres")
1106	{
1107	out[s]=list(toVdStrictSequence(C[s][1],d,v, Syzstring,s));
1108	}
1109	else
1110	{
1111	forout=toVdStrictSequence(C[s][1],d,v, Syzstring,s);
1112	list resolutionofA=forout[9];
1113	list resolutionofC=forout[10];
1114	forout=delete(forout,10);
1115	forout=delete(forout,9);
1116	out[s]=list(forout);
1117	for (i=1; i<=size(resolutionofC); i++)
1118	{
1119	out[s][1][5][i+1]=resolutionofC[i];//save the resolutions
1120	out[s][1][1][i+1]=resolutionofA[i];
1121	}
1122	}
1123	out[s][2]=list(list(out[s][1][3][1]));
1124	out[s][2][2]=list(unitmat(ncols(out[s][1][3][1])));
1125	out[s][2][3]=list(out[s][1][3][1]);
1126	out[s][2][4]=list(list());
1127	out[s][2][5]=list(list());
1128	out[s][2][6]=list(out[s][1][7][1]);
1129	out[s][2][7]=list(out[s][2][6][1]);
1130	out[s][2][8]=list(list());
1131	list resolutionofD;
1132	list resolutionofF;
1133	for (i=s-1; i>=2; i--)
1134	{
1135	C[i][2][5]=out[i+1][1][1][1];
1136	forout=toVdStrictSequences(C[i],d,out[i+1][1][6][1],Syzstring,s);
1137	if (Syzstring=="Sres")
1138	{
1139	resolutionofD=forout[3];//save the resolutions
1140	resolutionofF=forout[4];
1141	forout=delete(forout,4);
1142	forout=delete(forout,3);
1143	}
1144	out[i]=forout;
1145	if(Syzstring=="Sres")
1146	{
1147	for (j=2; j<=size(out[i+1][1][1]); j++)
1148	{
1149	out[i][2][5][j]=out[i+1][1][1][j];
1150	}
1151	for (j=1; j<=size(resolutionofD);j++)
1152	{
1153	out[i][1][1][j+1]=resolutionofD[j];
1154	out[i][1][5][j+1]=resolutionofF[j];
1155	}
1156	}
1157	}
1158	out[1]=list(list());//initalize our list
1159	C[1][2][5]=out[2][1][1][1];
1160	/Compute the last V_d-strict seqeunce/
1161	if (Syzstring=="Vdres")
1162	{
1163	out[1][2]=toVdStrictSequence(C[1][2],d,out[2][1][6][1],Syzstring,s,"J_Agiv");
1164	}
1165	else
1166	{
1167	forout=toVdStrictSequence(C[1][2],d,out[2][1][6][1],Syzstring,s,"J_Agiv");
1168	out[1][2]=delete(forout,9);
1169	list resolutionofA2=forout[9];
1170	for (i=1; i<=size(out[2][1][1]); i++)
1171	{
1172	/put the modules for the resolutions in the right spot/
1173	out[1][2][5][i]=out[2][1][1][i];
1174	}
1175	for (i=1; i<=size(resolutionofA2); i++)
1176	{
1177	out[1][2][1][i+1]=resolutionofA2[i];
1178	}
1179	}
1180	out[1][1][3]=list(out[1][2][1][1]);
1181	out[1][1][5]=list(out[1][2][1][1]);
1182	out[1][1][4]=list(unitmat(ncols(out[1][1][3][1])));
1183	out[1][1][7]=list(out[1][2][6][1]);
1184	out[1][1][8]=list(out[1][2][6][1]);
1185	out[1][1][1]=list(list());
1186	out[1][1][2]=list(list());
1187	out[1][1][6]=list(list());
1188	if (Syzstring=="Sres")
1189	{
1190	for (i=1; i<=size(out[1][2][1]); i++)
1191	{
1192	out[1][1][3][i]=out[1][2][1][i];
1193	out[1][1][5][i]=out[1][2][1][i];
1194	}
1195	}
1196	list Hi;
1197	for (i=1; i<=size(out); i++)
1198	{
1199	Hi[i]=list(out[i][1][5][1],out[i][1][8][1]);
1200	}
1201	list outall;
1202	outall[1]=out;
1203	outall[2]=Hi;
1204	return(outall);
1205	}
1206
1207	////////////////////////////////////////////////////////////////////////////////////
1208
1209	static proc toVdStrictSequence(list C,int n,intvec v,string Syzstring,int si,list #)
1210	{
1211	/this is the Algorithm 3.11 in [W2]/
1212	int omitemptylist;
1213	int lengthofres=si+n-1;
1214	int i,j,logens;
1215	def B=basering;
1216	matrix bi=slimgb(transpose(C[5]));
1217	/* Computation of a V_d-strict Groebner basis of C[5]:
1218	-if Syzstring=="Vdres" this is done using the method of weighted homogenization
1219	(Prop. 3.9 [OT])
1220	-else we use the homogenized Weyl algebra for Groebner basis computations
1221	(Prop 9.9 [OT]),
1222	in this case we already compute someresolutions (Thm. 9.10 [OT]) to omit
1223	extra Groebner basis computations later on*/
1224	int nr,nc;
1225	intvec i1,i2;
1226	if (Syzstring=="Vdres")
1227	{
1228	if(size(#)==0)
1229	{
1230	matrix J_C=VdStrictGB(C[5],n,list(v));
1231	}
1232	else
1233	{
1234	matrix J_C=C[5];//C[5] is already a V_d-strict Groebner basis
1235	}
1236	}
1237	else
1238	{
1239	if (size(#)==0)
1240	{
1241	matrix MC=C[5];
1242	def HomWeyl=makeHomogenizedWeyl(nvars(B) div 2, v);
1243	setring HomWeyl;
1244	matrix J_C=fetch(B,MC);
1245	J_C=nHomogenize(J_C);
1246	/*computation of V_d-strict resolution of C[5]->needed for proc
1247	VdstrictDoubleComplexes*/
1248	def ringofSyz=Sres(groebner(transpose(J_C)), lengthofres);
1249	setring ringofSyz;
1250	matrix J_C=transpose(matrix(RES[2]));
1251	J_C=subst(J_C,h,1);
1252	logens=ncols(J_C)+1;
1253	matrix zerom;
1254	list rofC;//will contain resolution of C
1255	for (i=3; i<=n+si+1; i++)
1256	{
1257	if (size(RES)>=i)
1258	{
1259	zerom=matrix(0,nrows(matrix(RES[i])),ncols(matrix(RES[i])));
1260	if (RES[i]!=zerom)
1261	{
1262	rofC[i-2]=(matrix(RES[i]));
1263
1264	if (i==3)
1265	{
1266	if (nrows(rofC[i-2])-logens+1!=nrows(J_C))
1267	{
1268	//build the resolution
1269	nr=nrows(J_C)+logens-1;
1270	nc=ncols(rofC[i-2]);
1271	rofC[i-2]=matrix(rofC[i-2],nr,nc);
1272	}
1273
1274	}
1275	if (i!=3)
1276	{
1277	if (nrows(rofC[i-2])-logens+1!=nrows(rofC[i-3]))
1278	{
1279	nr=nrows(rofC[i-3])+logens-1;
1280	nc=ncols(rofC[i-2]);
1281	rofC[i-2]=matrix(rofC[i-2],nr,nc);
1282	}
1283	}
1284	i1=intvec(logens..nrows(rofC[i-2]));
1285	i2=intvec(1..ncols(rofC[i-2]));
1286	rofC[i-2]=transpose(submat(rofC[i-2],i1,i2));
1287	logens=logens+ncols(rofC[i-2]);
1288	rofC[i-2]=subst(rofC[i-2],h,1);
1289	}
1290	else
1291	{
1292	rofC[i-2]=list();
1293	}
1294	}
1295	else
1296	{
1297	rofC[i-2]=list();
1298	}
1299	}
1300	if(size(rofC[1])==0)
1301	{
1302	omitemptylist=1;
1303	}
1304	setring B;
1305	matrix J_C=fetch(ringofSyz,J_C);
1306	if (omitemptylist!=1)
1307	{
1308	list rofC=fetch(ringofSyz,rofC);
1309	}
1310	omitemptylist=0;
1311	kill HomWeyl;
1312	kill ringofSyz;
1313	}
1314	else
1315	{
1316	matrix J_C=C[5];//C[5] is already a V_d-strict Groebner basis
1317	}
1318	}
1319	/* we compute a V_d-strict Groebner basis for C[3]*/
1320	matrix J_A=C[1];
1321	matrix f_CB=C[4];
1322	matrix f_ACB=transpose(concat(transpose(C[2]),transpose(f_CB)));
1323	matrix J_AC=divdr(f_ACB,C[3]);
1324	matrix P=matrixLift(J_AC * prodr(ncols(C[1]),ncols(C[5])) ,J_C);
1325	list storePi;
1326	matrix Pi[1][ncols(J_AC)];
1327	for (i=1; i<=nrows(J_C); i++)
1328	{
1329	for (j=1; j<=nrows(J_AC);j++)
1330	{
1331	Pi=Pi+P[i,j]*submat(J_AC,j,intvec(1..ncols(J_AC)));
1332	}
1333	storePi[i]=Pi;
1334	Pi=0;
1335	}
1336	/we compute the shift vector for C[1]/
1337	intvec m_a;
1338	list findMin;
1339	int comMin;
1340	for (i=1; i<=ncols(J_A); i++)
1341	{
1342	for (j=1; j<=size(storePi);j++)
1343	{
1344	if (storePi[j][1,i]!=0)
1345	{
1346	comMin=VdDeg(storePi[j]*prodr(ncols(J_A),ncols(C[5])),n,v);
1347	comMin=comMin-VdDeg(storePi[j][1,i],n,intvec(0));
1348	findMin[size(findMin)+1]=comMin;
1349	}
1350	}
1351	if (size(findMin)!=0)
1352	{
1353	m_a[i]=Min(findMin);
1354	findMin=list();
1355	}
1356	else
1357	{
1358	m_a[i]=0;
1359	}
1360	}
1361	matrix zero[ncols(J_A)][ncols(J_C)];
1362	matrix g_AB=concat(unitmat(ncols(J_A)),zero);
1363	matrix g_BC= transpose(concat(transpose(zero),transpose(unitmat(ncols(J_C)))));
1364	intvec m_b=m_a,v;
1365	/* computation of a V_d-strict Groebner basis of C[1] (and resolution if
1366	Syzstring=='Vdres') */
1367	if (Syzstring=="Vdres")
1368	{
1369	J_A=VdStrictGB(J_A,n,m_a);
1370	}
1371	else
1372	{
1373	def HomWeyl=makeHomogenizedWeyl(nvars(B) div 2, m_a);
1374	setring HomWeyl;
1375	matrix J_A=fetch(B,J_A);
1376	J_A=nHomogenize(J_A);
1377	def ringofSyz=Sres(groebner(transpose(J_A)),lengthofres);
1378	setring ringofSyz;
1379	matrix J_A=transpose(matrix(RES[2]));
1380	matrix zerom;
1381	J_A=subst(J_A,h,1);
1382	logens=ncols(J_A)+1;
1383	list rofA;
1384	for (i=3; i<=n+si+1; i++)
1385	{
1386	if (size(RES)>=i)
1387	{
1388	zerom=matrix(0,nrows(matrix(RES[i])),ncols(matrix(RES[i])));
1389	if (RES[i]!=zerom)
1390	{
1391	rofA[i-2]=matrix(RES[i]);// resolution for C[1]
1392	if (i==3)
1393	{
1394	if (nrows(rofA[i-2])-logens+1!=nrows(J_A))
1395	{
1396	nr=nrows(J_A)+logens-1;
1397	nc=ncols(rofA[i-2]);
1398	rofA[i-2]=matrix(rofA[i-2],nr,nc);
1399	}
1400	}
1401	if (i!=3)
1402	{
1403	if (nrows(rofA[i-2])-logens+1!=nrows(rofA[i-3]))
1404	{
1405	nr=nrows(rofA[i-3])+logens-1;
1406	nc=ncols(rofA[i-2]);
1407	rofA[i-2]=matrix(rofA[i-2],nr,nc);
1408	}
1409	}
1410	i1=intvec(logens..nrows(rofA[i-2]));
1411	i2=intvec(1..ncols(rofA[i-2]));
1412	rofA[i-2]=transpose(submat(rofA[i-2],i1,i2));
1413	logens=logens+ncols(rofA[i-2]);
1414	rofA[i-2]=subst(rofA[i-2],h,1);
1415	}
1416	else
1417	{
1418	rofA[i-2]=list();
1419	}
1420	}
1421	else
1422	{
1423	rofA[i-2]=list();
1424	}
1425	}
1426	if(size(rofA[1])==0)
1427	{
1428	omitemptylist=1;
1429	}
1430	setring B;
1431	J_A=fetch(ringofSyz,J_A);
1432	if (omitemptylist!=1)
1433	{
1434	list rofA=fetch(ringofSyz,rofA);
1435	}
1436	omitemptylist=0;
1437	kill HomWeyl;
1438	kill ringofSyz;
1439	}
1440	J_AC=transpose(storePi[1]);
1441	for (i=2; i<= size(storePi); i++)
1442	{
1443	J_AC=concat(J_AC, transpose(storePi[i]));
1444	}
1445	J_AC=transpose(concat(transpose(matrix(J_A,nrows(J_A),nrows(J_AC))),J_AC));
1446	list Vdstrict=(list(J_A),list(g_AB),list(J_AC),list(g_BC),list(J_C),list(m_a));
1447	Vdstrict[7]=list(m_b);
1448	Vdstrict[8]=list(v);
1449	if(Syzstring=="Sres")
1450	{
1451	Vdstrict[9]=rofA;
1452	if(size(#)==0)
1453	{
1454	Vdstrict[10]=rofC;
1455	}
1456	}
1457	return (Vdstrict);
1458	}
1459
1460	////////////////////////////////////////////////////////////////////////////////////
1461
1462	static proc toVdStrictSequences (list L,int d,intvec v,string Syzstring,int sizeL)
1463	{
1464	/* this is Argorithm 3.11 combined with Lemma 4.2 in [W2] for two short exact
1465	pieces.
1466	We asume that we are given two sequences of the form coker(L[i][1])->
1467	coker(L[i][3])->coker(L[i][5]) with differentials L[i][2] and L[i][4] such
1468	that L[1][3]=L[2][1].We are going to transform them to V_d-strict sequences
1469	J_D->J_A->J_F and J_A->J_B->J_C*/
1470	int omitemptylist;
1471	int lengthofres=sizeL+d-1;
1472	int logens;
1473	def B=basering;
1474	matrix J_F=L[1][5];
1475	matrix J_D=L[1][1];
1476	matrix f_FA=L[1][4];
1477	/*We find new presentations coker(J_DF) and coker(J_DFC) for L[1][4]=L[2][1]
1478	and L[2][4],resp. such that ncols(L[i][1])+ncols(L[i][5])=ncols(L[i][3]) */
1479	matrix f_DFA=transpose(concat(transpose(L[1][2]),transpose(f_FA)));
1480	matrix J_DF=divdr(f_DFA,L[1][3]);//coker(J_DF) is isomorphic to coker(L[2][1]);
1481	matrix J_C=L[2][5];
1482	matrix f_CB=L[2][4];
1483	matrix f_DFCB=transpose(concat(transpose(f_DFA*L[2][2]),transpose(f_CB)));
1484	matrix J_DFC=divdr(f_DFCB,L[2][3]);//coker(J_DFC) are coker(L[2][3)]) isomorphic
1485	/* find a shift vector on the range of J_F such that the first sequence is
1486	exact*/
1487	matrix P=matrixLift(J_DFC*prodr(ncols(J_DF),ncols(L[2][5])),J_C);
1488	list storePi;
1489	matrix Pi[1][ncols(J_DFC)];
1490	int i; int j;
1491	for (i=1; i<=nrows(J_C); i++)
1492	{
1493	for (j=1; j<=nrows(J_DFC);j++)
1494	{
1495	Pi=Pi+P[i,j]*submat(J_DFC,j,intvec(1..ncols(J_DFC)));
1496	}
1497	storePi[i]=Pi;
1498	Pi=0;
1499	}
1500	intvec m_a;
1501	list findMin;
1502	list noMin;
1503	int comMin;
1504	int nr,nc;
1505	intvec i1,i2;
1506	for (i=1; i<=ncols(J_DF); i++)
1507	{
1508	for (j=1; j<=size(storePi);j++)
1509	{
1510	if (storePi[j][1,i]!=0)
1511	{
1512	comMin=VdDeg(storePi[j]*prodr(ncols(J_DF),ncols(J_C)),d,v);
1513	comMin=comMin-VdDeg(storePi[j][1,i],d,intvec(0));
1514	findMin[size(findMin)+1]=comMin;
1515	}
1516	}
1517	if (size(findMin)!=0)
1518	{
1519	m_a[i]=Min(findMin);// shift vector for L[2][1]
1520	findMin=list();
1521	noMin[i]=0;
1522	}
1523	else
1524	{
1525	noMin[i]=1;
1526	}
1527	}
1528	if (size(m_a) < ncols(J_DF))
1529	{
1530	m_a[ncols(J_DF)]=0;
1531	}
1532	intvec m_f=m_a[ncols(J_D)+1..size(m_a)];
1533	/* Computation of a V_d-strict Groebner basis of J_F=L[1][5]:
1534	if Syzstring=="Vdres" this is done using the method of weighted homogenization
1535	(Prop. 3.9 [OT])
1536	else we use the homogenized Weyl algerba for Groebner basis computations
1537	(Prop 9.9 [OT]), in this case we already compute resolutions
1538	(Thm. 9.10 in [OT]) to omit extra Groebner basis computations later on*/
1539	if (Syzstring=="Vdres")
1540	{
1541	J_F=VdStrictGB(J_F,d,m_f);
1542	}
1543	else
1544	{
1545	def HomWeyl=makeHomogenizedWeyl(nvars(B) div 2, m_f);
1546	setring HomWeyl;
1547	matrix J_F=fetch(B,J_F);
1548	J_F=nHomogenize(J_F);
1549	def ringofSyz=Sres(groebner(transpose(J_F)),lengthofres);
1550	setring ringofSyz;
1551	matrix J_F=transpose(matrix(RES[2]));
1552	J_F=subst(J_F,h,1);
1553	logens=ncols(J_F)+1;
1554	list rofF;
1555	for (i=3; i<=d+sizeL+1; i++)
1556	{
1557	if (size(RES)>=i)
1558	{
1559	if (RES[i]!=matrix(0,nrows(matrix(RES[i])),ncols(matrix(RES[i]))))
1560	{
1561	rofF[i-2]=(matrix(RES[i]));// resolution for J_F
1562	if (i==3)
1563	{
1564	if (nrows(rofF[i-2])-logens+1!=nrows(J_F))
1565	{
1566	nr=nrows(J_F)+logens-1;
1567	nc=ncols(rofF[i-2]);
1568	rofF[i-2]=matrix(rofF[i-2],nr,nc);
1569	}
1570	}
1571	if (i!=3)
1572	{
1573	if (nrows(rofF[i-2])-logens+1!=nrows(rofF[i-3]))
1574	{
1575	nr=nrows(rofF[i-3])+logens-1;
1576	rofF[i-2]=matrix(rofF[i-2],nr,ncols(rofF[i-2]));
1577	}
1578	}
1579	i1=intvec(logens..nrows(rofF[i-2]));
1580	i2=intvec(1..ncols(rofF[i-2]));
1581	rofF[i-2]=transpose(submat(rofF[i-2],i1,i2));
1582	logens=logens+ncols(rofF[i-2]);
1583	rofF[i-2]=subst(rofF[i-2],h,1);
1584	}
1585	else
1586	{
1587	rofF[i-2]=list();
1588	}
1589	}
1590	else
1591	{
1592	rofF[i-2]=list();
1593	}
1594	}
1595	if(size(rofF[1])==0)
1596	{
1597	omitemptylist=1;
1598	}
1599	setring B;
1600	J_F=fetch(ringofSyz,J_F);
1601	if (omitemptylist!=1)
1602	{
1603	list rofF=fetch(ringofSyz,rofF);
1604	}
1605	omitemptylist=0;
1606	kill HomWeyl;
1607	kill ringofSyz;
1608	}
1609	/find shift vectors on the range of J_D/
1610	P=matrixLift(J_DF * prodr(ncols(L[1][1]),ncols(L[1][5])) ,J_F);
1611	list storePinew;
1612	matrix Pidf[1][ncols(J_DF)];
1613	for (i=1; i<=nrows(J_F); i++)
1614	{
1615	for (j=1; j<=nrows(J_DF);j++)
1616	{
1617	Pidf=Pidf+P[i,j]*submat(J_DF,j,intvec(1..ncols(J_DF)));
1618	}
1619	storePinew[i]=Pidf;
1620	Pidf=0;
1621	}
1622	intvec m_d;
1623	for (i=1; i<=ncols(J_D); i++)
1624	{
1625	for (j=1; j<=size(storePinew);j++)
1626	{
1627	if (storePinew[j][1,i]!=0)
1628	{
1629	comMin=VdDeg(storePinew[j]*prodr(ncols(J_D),ncols(L[1][5])),d,m_f);
1630	comMin=comMin-VdDeg(storePinew[j][1,i],d,intvec(0));
1631	findMin[size(findMin)+1]=comMin;
1632	}
1633	}
1634	if (size(findMin)!=0)
1635	{
1636	if (noMin[i]==0)
1637	{
1638	m_d[i]=Min(insert(findMin,m_a[i]));
1639	m_a[i]=m_d[i];
1640	}
1641	else
1642	{
1643	m_d[i]=Min(findMin);
1644	m_a[i]=m_d[i];
1645	}
1646	}
1647	else
1648	{
1649	m_d[i]=m_a[i];
1650	}
1651	findMin=list();
1652	}
1653	/* compute a V_d-strict Groebner basis (and resolution of J_D if
1654	Syzstring!='Vdres') for J_D*/
1655	if (Syzstring=="Vdres")
1656	{
1657	J_D=VdStrictGB(J_D,d,m_d);
1658	}
1659	else
1660	{
1661	def HomWeyl=makeHomogenizedWeyl(nvars(B) div 2, m_d);
1662	setring HomWeyl;
1663	matrix J_D=fetch(B,J_D);
1664	J_D=nHomogenize(J_D);
1665	def ringofSyz=Sres(groebner(transpose(J_D)),lengthofres);
1666	setring ringofSyz;
1667	matrix J_D=transpose(matrix(RES[2]));
1668	J_D=subst(J_D,h,1);
1669	logens=ncols(J_D)+1;
1670	list rofD;
1671	for (i=3; i<=d+sizeL+1; i++)
1672	{
1673	if (size(RES)>=i)
1674	{
1675	if (RES[i]!=matrix(0,nrows(matrix(RES[i])),ncols(matrix(RES[i]))))
1676	{
1677	rofD[i-2]=(matrix(RES[i]));// resolution for J_D
1678	if (i==3)
1679	{
1680	if (nrows(rofD[i-2])-logens+1!=nrows(J_D))
1681	{
1682	nr=nrows(J_D)+logens-1;
1683	rofD[i-2]=matrix(rofD[i-2],nr,ncols(rofD[i-2]));
1684	}
1685	}
1686	if (i!=3)
1687	{
1688	if (nrows(rofD[i-2])-logens+1!=nrows(rofD[i-3]))
1689	{
1690	nr=nrows(rofD[i-3])+logens-1;
1691	rofD[i-2]=matrix(rofD[i-2],nr,ncols(rofD[i-2]));
1692	}
1693	}
1694	i1=intvec(logens..nrows(rofD[i-2]));
1695	i2=intvec(1..ncols(rofD[i-2]));
1696	rofD[i-2]=transpose(submat(rofD[i-2],i1,i2));
1697	logens=logens+ncols(rofD[i-2]);
1698	rofD[i-2]=subst(rofD[i-2],h,1);
1699	}
1700	else
1701	{
1702	rofD[i-2]=list();
1703	}
1704	}
1705	else
1706	{
1707	rofD[i-2]=list();
1708	}
1709	}
1710	if(size(rofD[1])==0)
1711	{
1712	omitemptylist=1;
1713	}
1714	setring B;
1715	J_D=fetch(ringofSyz,J_D);
1716	if (omitemptylist!=1)
1717	{
1718	list rofD=fetch(ringofSyz,rofD);
1719	}
1720	omitemptylist=0;
1721	kill HomWeyl;
1722	kill ringofSyz;
1723	}
1724	/* compute new matrices for J_A and J_B such that their rows form a V_d-strict
1725	Groebner basis and nrows(J_A)=nrows(J_D)+nrows(J_F) and
1726	nrows(J_B)=nrows(J_A)+nrows(J_C)*/
1727	J_DF=transpose(storePinew[1]);
1728	for (i=2; i<=nrows(J_F); i++)
1729	{
1730	J_DF=concat(J_DF,transpose(storePinew[i]));
1731	}
1732	J_DF=transpose(concat(transpose(matrix(J_D,nrows(J_D),nrows(J_DF))),J_DF));
1733	J_DFC=transpose(storePi[1]);
1734	for (i=2; i<=nrows(J_C); i++)
1735	{
1736	J_DFC=concat(J_DFC,transpose(storePi[i]));
1737	}
1738	J_DFC=transpose(concat(transpose(matrix(J_DF,nrows(J_DF),nrows(J_DFC))),J_DFC));
1739	intvec m_b=m_a,v;
1740	matrix zero[ncols(J_D)][ncols(J_F)];
1741	matrix g_DA=concat(unitmat(ncols(J_D)),zero);
1742	matrix g_AF=transpose(concat(transpose(zero),unitmat(ncols(J_F))));
1743	matrix zero1[ncols(J_DF)][ncols(J_C)];
1744	matrix g_AB=concat(unitmat(ncols(J_DF)),zero1);
1745	matrix g_BC=transpose(concat(transpose(zero1),unitmat(ncols(J_C))));
1746	list out;
1747	out[1]=list(list(J_D),list(g_DA),list(J_DF),list(g_AF),list(J_F));
1748	out[1]=out[1]+list(list(m_d),list(m_a),list(m_f));
1749	out[2]=list(list(J_DF),list(g_AB),list(J_DFC),list(g_BC),list(J_C));
1750	out[2]=out[2]+list(list(m_a),list(m_b),list(v));
1751	if (Syzstring=="Sres")
1752	{
1753	out[3]=rofD;
1754	out[4]=rofF;
1755	}
1756	return(out);
1757	}
1758
1759	////////////////////////////////////////////////////////////////////////////////////
1760
1761	static proc VdStrictDoubleComplexes(list L,int d,string Syzstring)
1762	{
1763	/* We compute V_d-strict resolutions over the V_d-strict short exact pieces from
1764	the procedure shortExactPiecesToVdStrict.
1765	We use Algorithms 3.14 and 3.15 in [W2]*/
1766	int i,k,c,j,l,totaldeg,comparedegs,SBcom,verk;
1767	intvec fordegs;
1768	intvec n_b,i1,i2;
1769	matrix rem,forML,subm,zerom,unitm,subm2;
1770	matrix J_B;
1771	list store;
1772	int t=size(L)+d;
1773	int vd1,vd2,nr,nc;
1774	def B=basering;
1775	int n=nvars(B) div 2;
1776	intvec v;
1777	list forhW;
1778	if (Syzstring=="Sres")
1779	{
1780	/*we already computed some of the resolutions in the procedure
1781	shortExactPiecesToVdStrict*/
1782	matrix Pold,Pnew,Picombined; intvec containsndeg; matrix Pinew;
1783	for (k=1; k<=(size(L)+d-1); k++)
1784	{
1785	L[1][1][1][k+1]=list();
1786	L[1][1][2][k+1]=list();
1787	L[1][1][6][k+1]=list();
1788	}
1789	L[1][1][6][size(L)+d+1]=list();
1790	matrix mem;
1791	for (i=2; i<=d+size(L)+1; i++)
1792	{;
1793	v=0;
1794	if(size(L[1][1][3][i-1])!=0)
1795	{
1796	if(i!=d+size(L)+1)
1797	{
1798	/horizontal differential/
1799	L[1][1][4][i-1]=unitmat(nrows(L[1][1][3][i-1]));
1800	}
1801	for (j=1; j<=nrows(L[1][1][3][i-1]); j++)
1802	{
1803	mem=submat(L[1][1][3][i-1],j,intvec(1..ncols(L[1][1][3][i-1])));
1804	v[j]=VdDeg(mem,d,L[1][1][7][i-1]);
1805	}
1806	L[1][1][7][i]=v;//new shift vector
1807	L[1][1][8][i]=v;
1808	L[1][2][6][i]=v;
1809	}
1810	else
1811	{
1812	if (i!=d+size(L)+1)
1813	{
1814	L[1][1][4][i-1]=list();
1815	}
1816	L[1][1][7][i]=list();
1817	L[1][1][8][i]=list();
1818	L[1][2][6][i]=list();
1819	}
1820	}
1821	if (size(L[1][1][3][d+size(L)])!=0)
1822	{
1823	/horizontal differential/
1824	L[1][1][4][d+size(L)]=unitmat(nrows(L[1][1][3][d+size(L)]));
1825	}
1826	else
1827	{
1828	L[1][1][4][d+size(L)]=list();
1829	}
1830	for (k=1; k<size(L); k++)
1831	{
1832	/* We build a V_d-strict resolution for coker(L[k][2][1][1])->
1833	coker(L[k][2][3][1])->coker(L[k][2][5][1]) using the resolution
1834	obtained for coker(L[k][1][3][1]).
1835	L[k][2][i][j] will be the jth module in the resolution of L[k][2][i][1]
1836	for i=1,3,5.
1837	L[k][2][i+5][j] will be the jth shift vector in the resolution of
1838	L[k][2][i][1](this holds also for the case Syzstring=="Vdres")*/
1839	for (i=2; i<=d+size(L); i++)
1840	{
1841	v=0;
1842	if (size(L[k][2][5][i-1])!=0)
1843	{
1844	for (j=1; j<=nrows(L[k][2][5][i-1]); j++)
1845	{
1846	i1=intvec(1..ncols(L[k][2][5][i-1]));
1847	mem=submat(L[k][2][5][i-1],j,i1);
1848	v[j]=VdDeg(mem,d,L[k][2][8][i-1]);
1849	}
1850	/next shift vector in th resolution of coker(L[k][2][5][1])/
1851	L[k][2][8][i]=v;
1852	}
1853	else
1854	{
1855	L[k][2][8][i]=list();
1856	}
1857	/* we build step by step a resolution for coker(L[k][2][5][1]) using
1858	the resolutions of coker(L[k][2][1][1]) and coker(L[k][2][5][1])*/
1859	if (size(L[k][2][5][i])!=0)
1860	{
1861	if (size(L[k][2][1][i])!=0 or size(L[k][2][1][i-1])!=0)
1862	{
1863	L[k][2][3][i]=transpose(syz(transpose(L[k][2][3][i-1])));
1864	nr= nrows(L[k][2][1][i-1]);
1865	nc=ncols(L[k][2][5][i]);
1866	Pold=matrixLift(L[k][2][3][i]*prodr(nr,nc), L[k][2][5][i]);
1867	matrix Pi[1][ncols(L[k][2][3][i])];
1868	for (l=1; l<=nrows(L[k][2][5][i]); l++)
1869	{
1870	for (j=1; j<=nrows(L[k][2][3][i]); j++)
1871	{
1872	i2=intvec(1..ncols(L[k][2][3][i]));
1873	Pi=Pi+Pold[l,j]*submat(L[k][2][3][i],j,i2);
1874	}
1875	if (l==1)
1876	{
1877	Picombined=transpose(Pi);
1878	}
1879	else
1880	{
1881	Picombined=concat(Picombined,transpose(Pi));
1882	}
1883	Pi=0;
1884	}
1885	kill Pi;
1886	Picombined=transpose(Picombined);
1887	if (size(L[k][2][1][i])!=0)
1888	{
1889	if (i==2)
1890	{
1891	containsndeg=(0:ncols(L[k][2][1][1]));
1892	}
1893	containsndeg=nDeg(L[k][2][1][i-1],containsndeg);
1894	forhW=list(L[k][2][6][i],containsndeg);
1895	def HomWeyl=makeHomogenizedWeyl(n,forhW);
1896	setring HomWeyl;
1897	list L=fetch(B,L);
1898	matrix M=L[k][2][1][i];
1899	list forM=nHomogenize(M,containsndeg,1);
1900	M=forM[1];
1901	totaldeg=forM[2];
1902	kill forM;
1903	matrix Maorig=fetch(B,Picombined);
1904	matrix Ma=submat(Maorig,(1..nrows(Maorig)),(1..ncols(M)));
1905	matrix mem,subm,zerom;
1906	matrix Pinew;
1907	M=transpose(M);
1908	SBcom=0;
1909	for (l=1; l<=nrows(Ma); l++)
1910	{
1911	zerom=matrix(0,1,(ncols(Maorig)-ncols(Ma)));
1912	i1=(ncols(Ma)+1..ncols(Maorig));
1913	if (submat(Maorig,l,i1)==zerom)
1914	{
1915	for (cc=1; cc<=ncols(Ma); cc++)
1916	{
1917	Maorig[l,cc]=0;
1918	}
1919	}
1920	i2=(ncols(Ma)+1..ncols(Maorig));
1921	i1=(1..ncols(Ma));
1922	if (VdDeg(submat(Maorig,l,i1),d,L[k][2][6][i])>
1923	VdDeg(submat(Maorig,l,i2),d,L[k][2][8][i]) and
1924	submat(Maorig,l,i1)!=matrix(0,1,ncols(Ma)))
1925	{
1926	/*V_d-Grad is to big--> we make it smaller using
1927	Vdnormal form computations*/
1928	if (SBcom==0)
1929	{
1930	M=slimgb(M);
1931	SBcom=1;
1932	}
1933	//print("Reduzierung des V_d-Grades(Stelle1)");
1934	i2=(ncols(Ma)+1..ncols(Maorig));
1935	vd1=VdDeg(submat(Maorig,l,i2),d,L[k][2][8][i]);
1936	mem=submat(Ma,l,(1..ncols(Ma)));
1937	mem=nHomogenize(mem,containsndeg);
1938	mem=h^totaldeg*mem;
1939	mem=transpose(mem);
1940	mem=reduce(mem,groebner(M));
1941	matrix jt=transpose(subst(mem,h,1));
1942	setring B;
1943	matrix jt=fetch(HomWeyl,jt);
1944	matrix need=fetch(HomWeyl,Maorig);
1945	need=submat(need,l,(1..ncols(need)));
1946	i1=L[k][2][6][i];
1947	i2=L[k][2][8][i];
1948	jt=VdNormalForm(need,L[k][2][1][i],d,i1,i2);
1949	setring HomWeyl;
1950	mem=fetch(B,jt);
1951	mem=transpose(mem);
1952	if (l==1)
1953	{
1954	Pinew=mem;
1955	}
1956	else
1957	{
1958	Pinew=concat(Pinew,mem);
1959	}
1960	vd2=VdDeg(transpose(mem),d,L[k][2][6][i]);
1961	if (vd2>vd1 and mem!=matrix(0,nrows(mem),ncols(mem)))
1962	{//should not happen!!
1963	//print("Reduzierung fehlgeschlagen!!(Stelle1)");
1964	}
1965	}
1966	else
1967	{
1968	if (l==1)
1969	{
1970	Pinew=transpose(submat(Ma,l,(1..ncols(Ma))));
1971	}
1972	else
1973	{
1974	subm=transpose(submat(Ma,l,(1..ncols(Ma))));
1975	Pinew=concat(Pinew,subm);
1976	}
1977	}
1978	}
1979	Pinew=subst(Pinew,h,1);
1980	Pinew=transpose(Pinew);
1981	setring B;
1982	Pinew=fetch(HomWeyl,Pinew);
1983	kill HomWeyl;
1984	L[k][2][3][i]=concat(Pinew,L[k][2][5][i]);
1985	subm=transpose(L[k][2][3][i]);
1986	subm=concat(transpose(L[k][2][1][i]),subm);
1987	L[k][2][3][i]=transpose(subm);
1988	}
1989	else
1990	{
1991	L[k][2][3][i]=Picombined;
1992	}
1993	L[k+1][1][1][i]=L[k][2][5][i];
1994	nr=nrows(L[k][2][1][i-1]);
1995	nc=ncols(L[k][2][5][i]);
1996	L[k][2][2][i]=concat(unitmat(nr),matrix(0,nr,nc));
1997	L[k][2][4][i]=prodr(nrows(L[k][2][1][i-1]),nc);
1998	v=L[k][2][6][i],L[k][2][8][i];
1999	L[k][2][7][i]=v;
2000	L[k+1][1][6][i]=L[k][2][8][i];
2001	}
2002	else
2003	{
2004	L[k][2][3][i]=L[k][2][5][i];
2005	L[k][2][2][i]=list();
2006	L[k][2][7][i]=L[k][2][8][i];
2007	L[k][2][4][i]=unitmat(nrows(L[k][2][5][i-1]));
2008	L[k+1][1][6][i]=L[k][2][8][i];
2009	L[k+1][1][1][i]=L[k][2][5][i];
2010	}
2011	}
2012	else
2013	{
2014	if (size(L[k][2][1][i])!=0)
2015	{
2016	if (size(L[k][2][5][i-1])!=0)
2017	{
2018	nr=nrows(L[k][2][5][i-1]);
2019	L[k][2][3][i]=concat(L[k][2][1][i],matrix(0,1,nr));
2020	v=L[k][2][6][i],L[k][2][8][i];
2021	L[k][2][7][i]=v;
2022	nc=nrows(L[k][2][1][i-1]);
2023	L[k][2][2][i]=concat(unitmat(nc),matrix(0,nc,nr));
2024	L[k][2][4][i]=prodr(nrows(L[k][2][1][i-1]),nr);
2025	}
2026	else
2027	{
2028	L[k][2][3][i]=L[k][2][1][i];
2029	L[k][2][7][i]=L[k][2][6][i];
2030	L[k][2][2][i]=unitmat(nrows(L[k][2][1][i-1]));
2031	L[k][2][4][i]=list();
2032	}
2033	L[k+1][1][1][i]=L[k][2][5][i];
2034	L[k+1][1][6][i]=L[k][2][8][i];
2035	}
2036	else
2037	{
2038	L[k][2][3][i]=list();
2039	if (size(L[k][2][6][i])!=0)
2040	{
2041	if (size(L[k][2][8][i])!=0)
2042	{
2043	v=L[k][2][6][i],L[k][2][8][i];
2044	L[k][2][7][i]=v;
2045	nr=nrows(L[k][2][1][i-1]);
2046	nc=nrows(L[k][2][5][i-1]);
2047	L[k][2][2][i]=concat(unitmat(nc),matrix(0,nr,nc));
2048	L[k][2][4][i]=prodr(nr,nrows(L[k][2][5][i-1]));
2049	}
2050	else
2051	{
2052	L[k][2][7][i]=L[k][2][6][i];
2053	L[k][2][2][i]=unitmat(nrows(L[k][2][1][i-1]));
2054	L[k][2][4][i]=list();
2055	}
2056	}
2057	else
2058	{
2059	if (size(L[k][2][8][i])!=0)
2060	{
2061	L[k][2][7][i]=L[k][2][8][i];
2062	L[k][2][2][i]=list();
2063	L[k][2][4][i]=unitmat(nrows(L[k][2][5][i-1]));
2064	}
2065	else
2066	{
2067	L[k][2][7][i]=list();
2068	L[k][2][2][i]=list();
2069	L[k][2][4][i]=list();
2070	}
2071	}
2072	L[k+1][1][1][i]=L[k][2][5][i];
2073	L[k+1][1][6][i]=L[k][2][8][i];
2074	}
2075	}
2076	}
2077	i=d+size(L)+1;
2078	v=0;
2079	if (size(L[k][2][5][i-1])!=0)
2080	{
2081	for (j=1; j<=nrows(L[k][2][5][i-1]); j++)
2082	{
2083	mem=submat(L[k][2][5][i-1],j,intvec(1..ncols(L[k][2][5][i-1])));
2084	v[j]=VdDeg(mem,d,L[k][2][8][i-1]);
2085	}
2086	L[k][2][8][i]=v;
2087	if (size(L[k][2][6][i])!=0)
2088	{
2089	v=L[k][2][6][i],L[k][2][8][i];
2090	L[k][2][7][i]=v;
2091	}
2092	else
2093	{
2094	L[k][2][7][i]=L[k][2][8][i];
2095	}
2096	}
2097	else
2098	{
2099	L[k][2][8][i]=list();
2100	L[k][2][7][i]=L[k][2][6][i];
2101	}
2102	L[k+1][1][6][i]=L[k][2][8][i];
2103	/* now we build V_d-strict resolutions for the sequences
2104	coker(L[k+1][1][1][1])->coker(L[k+1][1][3][1])->coker(L[k+1][1][5][i])
2105	using the resolutions for coker(L[k][2][5][1]) we just obtained
2106	(works exactly the same as above)*/
2107	for (i=2; i<=d+size(L); i++)
2108	{
2109	v=0;
2110	if (size(L[k+1][1][5][i-1])!=0)
2111	{
2112	for (j=1; j<=nrows(L[k+1][1][5][i-1]); j++)
2113	{
2114	i1=intvec(1..ncols(L[k+1][1][5][i-1]));
2115	mem=submat(L[k+1][1][5][i-1],j,i1);
2116	v[j]=VdDeg(mem,d,L[k+1][1][8][i-1]);
2117	}
2118	L[k+1][1][8][i]=v;
2119	}
2120	else
2121	{
2122	L[k+1][1][8][i]=list();
2123	}
2124	if (size(L[k+1][1][5][i])!=0)
2125	{
2126	if (size(L[k+1][1][1][i])!=0 or size(L[k+1][1][1][i-1])!=0)
2127	{
2128	L[k+1][1][3][i]=transpose(syz(transpose(L[k+1][1][3][i-1])));
2129	nr=nrows(L[k+1][1][1][i-1]);
2130	nc=ncols(L[k+1][1][5][i]);
2131	Pold=matrixLift(L[k+1][1][3][i]*prodr(nr,nc),L[k+1][1][5][i]);
2132	matrix Pi[1][ncols(L[k+1][1][3][i])];
2133	for (l=1; l<=nrows(L[k+1][1][5][i]); l++)
2134	{
2135	for (j=1; j<=nrows(L[k+1][1][3][i]); j++)
2136	{
2137	i2=intvec(1..ncols(L[k+1][1][3][i]));
2138	Pi=Pi+Pold[l,j]*submat(L[k+1][1][3][i],j,i2);
2139	}
2140	if (l==1)
2141	{
2142	Picombined=transpose(Pi);
2143	}
2144	else
2145	{
2146	Picombined=concat(Picombined,transpose(Pi));
2147	}
2148	Pi=0;
2149	}
2150	kill Pi;
2151	Picombined=transpose(Picombined);
2152	if(size(L[k+1][1][1][i])!=0)
2153	{
2154	if (i==2)
2155	{
2156	containsndeg=(0:ncols(L[k+1][1][1][i-1]));
2157	}
2158	containsndeg=nDeg(L[k+1][1][1][i-1],containsndeg);
2159	forhW=list(L[k+1][1][6][i], containsndeg);
2160	def HomWeyl=makeHomogenizedWeyl(n,forhW);
2161	setring HomWeyl;
2162	list L=fetch(B,L);
2163	matrix M=L[k+1][1][1][i];
2164	list forM=nHomogenize(M,containsndeg,1);
2165	M=forM[1];
2166	totaldeg=forM[2];
2167	kill forM;
2168	matrix Maorig=fetch(B,Picombined);
2169	matrix Ma=submat(Maorig,(1..nrows(Maorig)),(1..ncols(M)));
2170	Ma=nHomogenize(Ma,containsndeg);
2171	matrix mem,subm,zerom,subm2;
2172	matrix Pinew;
2173	M=transpose(M);
2174	SBcom=0;
2175	for (l=1; l<=nrows(Ma); l++)
2176	{
2177	i2=(ncols(Ma)+1..ncols(Maorig));
2178	nc=ncols(Maorig)-ncols(Ma);
2179	if (submat(Maorig,l,i2)==matrix(0,1,nc))
2180	{
2181	for (cc=1; cc<=ncols(Ma); cc++)
2182	{
2183	Maorig[l,cc]=0;
2184	}
2185	}
2186	i1=(1..ncols(Ma));
2187	i2=L[k+1][1][8][i];
2188	subm=submat(Maorig,l,i1);
2189	subm2=submat(Maorig,l,(ncols(Ma)+1..ncols(Maorig)));
2190	if (VdDeg(subm,d,L[k+1][1][6][i])>VdDeg(subm2,d,i2)
2191	and subm!=matrix(0,1,ncols(Ma)))
2192	{
2193	//print("Reduzierung des Vd-Grades (Stelle2)");
2194	if (SBcom==0)
2195	{
2196	M=slimgb(M);
2197	SBcom=1;
2198	}
2199	vd1=VdDeg(subm2,d,L[k+1][1][8][i]);
2200	mem=submat(Ma,l,(1..ncols(Ma)));
2201	mem=nHomogenize(mem,containsndeg);
2202	mem=h^totaldeg*mem;
2203	mem=transpose(mem);
2204	mem=reduce(mem, groebner(M));
2205	if (l==1)
2206	{
2207	Pinew=mem;
2208	}
2209	else
2210	{
2211	Pinew=concat(Pinew,mem);
2212	}
2213	vd2=VdDeg(transpose(mem),d,L[k+1][1][6][i]);
2214	if (vd2>vd1 and mem!=matrix(0,nrows(mem),ncols(mem)))
2215	{//should not happen
2216	//print("Reduzierung fehlgeschlagen!!!!(Stelle2)");
2217	}
2218	}
2219	else
2220	{
2221	if (l==1)
2222	{
2223	Pinew=transpose(submat(Ma,l,(1..ncols(Ma))));
2224	}
2225	else
2226	{
2227	subm=transpose(submat(Ma,l,(1..ncols(Ma))));
2228	Pinew=concat(Pinew,subm);
2229	}
2230	}
2231	}
2232	Pinew=subst(Pinew,h,1);
2233	Pinew=transpose(Pinew);
2234	setring B;
2235	Pinew=fetch(HomWeyl,Pinew);
2236	kill HomWeyl;
2237	L[k+1][1][3][i]=concat(Pinew,L[k+1][1][5][i]);
2238	subm=transpose(L[k+1][1][1][i]);
2239	subm2=transpose(L[k+1][1][3][i]);
2240	L[k+1][1][3][i]=transpose(concat(subm,subm2));
2241	}
2242	else
2243	{
2244	L[k+1][1][3][i]=Picombined;
2245	}
2246	L[k+1][2][1][i]=L[k+1][1][3][i];
2247	nr=nrows(L[k+1][1][1][i-1]);
2248	nc=ncols(L[k+1][1][5][i]);
2249	L[k+1][1][2][i]=concat(unitmat(nr),matrix(0,nr,nc));
2250	L[k+1][1][4][i]=prodr(nr,nc);
2251	v=L[k+1][1][6][i],L[k+1][1][8][i];
2252	L[k+1][1][7][i]=v;
2253	L[k+1][2][6][i]=L[k+1][1][7][i];
2254	}
2255	else
2256	{
2257	L[k+1][1][3][i]=L[k+1][1][5][i];
2258	L[k+1][1][2][i]=list();
2259	L[k+1][1][4][i]=unitmat(nrows(L[k+1][1][5][i-1]));
2260	L[k+1][1][7][i]=L[k+1][1][8][i];
2261	L[k+1][2][6][i]=L[k+1][1][7][i];
2262	L[k+1][2][1][i]=L[k+1][1][3][i];
2263	}
2264	}
2265	else
2266	{
2267	if (size(L[k+1][1][1][i])!=0)
2268	{
2269	if (size(L[k+1][1][5][i-1])!=0)
2270	{
2271	zerom=matrix(0,1,nrows(L[k+1][1][5][i-1]));
2272	L[k+1][1][3][i]=concat(L[k+1][1][1][i],zerom);
2273	v=L[k+1][1][6][i],L[k+1][1][8][i];
2274	L[k+1][1][7][i]=v;
2275	nr=nrows(L[k+1][1][1][i-1]);
2276	nc=nrows(L[k+1][1][5][i-1]);
2277	L[k+1][1][2][i]=concat(unitmat(nr),matrix(0,nr,nc));
2278	L[k+1][1][4][i]=prodr(nr,nc);
2279	}
2280	else
2281	{
2282	L[k+1][1][3][i]=L[k+1][1][1][i];
2283	L[k+1][1][7][i]=L[k+1][1][6][i];
2284	L[k+1][1][2][i]=unitmat(nrows(L[k+1][1][1][i-1]));
2285	L[k+1][1][4][i]=list();
2286	}
2287	L[k+1][2][1][i]=L[k+1][1][3][i];
2288	L[k+1][2][6][i]=L[k+1][1][7][i];
2289	}
2290	else
2291	{
2292	L[k+1][1][3][i]=list();
2293	if (size(L[k+1][1][6][i])!=0)
2294	{
2295	if (size(L[k+1][1][8][i])!=0)
2296	{
2297	v=L[k+1][1][6][i],L[k+1][1][8][i];
2298	L[k+1][1][7][i]=v;
2299	nr=nrows(L[k+1][1][1][i-1]);
2300	nc=nrows(L[k+1][1][5][i-1]);
2301	L[k+1][1][2][i]=concat(unitmat(nr),matrix(0,nr,nc));
2302	L[k+1][1][4][i]=prodr(nr,nrows(L[k+1][1][5][i-1]));
2303	}
2304	else
2305	{
2306	L[k+1][1][7][i]=L[k+1][1][6][i];
2307	L[k+1][1][2][i]=unitmat(nrows(L[k+1][1][1][i-1]));
2308	L[k+1][1][4][i]=list();
2309	}
2310	}
2311	else
2312	{
2313	if (size(L[k+1][1][8][i])!=0)
2314	{
2315	L[k+1][1][7][i]=L[k+1][1][8][i];
2316	L[k+1][1][2][i]=list();
2317	L[k+1][1][4][i]=unitmat(nrows(L[k+1][1][5][i-1]));
2318	}
2319	else
2320	{
2321	L[k+1][1][7][i]=list();
2322	L[k+1][1][2][i]=list();
2323	L[k+1][1][4][i]=list();
2324	}
2325	}
2326
2327	L[k+1][2][1][i]=L[k+1][1][3][i];
2328	L[k+1][2][6][i]=L[k+1][1][7][i];
2329	}
2330	}
2331	}
2332	i=size(L)+d+1;
2333	v=0;
2334	if (size(L[k+1][1][5][i-1])!=0)
2335	{
2336	for (j=1; j<=nrows(L[k+1][1][5][i-1]); j++)
2337	{
2338	i1=intvec(1..ncols(L[k+1][1][5][i-1]));
2339	mem=submat(L[k+1][1][5][i-1],j,i1);
2340	v[j]=VdDeg(mem,d,L[k+1][1][8][i-1]);
2341	}
2342	L[k+1][1][8][i]=v;
2343	if (size(L[k+1][1][6][i])!=0)
2344	{
2345	v=L[k+1][1][6][i],L[k+1][1][8][i];
2346	L[k+1][1][7][i]=v;
2347	}
2348	else
2349	{
2350	L[k+1][1][7][i]=L[k+1][1][8][i];
2351	}
2352	}
2353	else
2354	{
2355	L[k+1][1][8][i]=list();
2356	L[k+1][1][7][i]=L[k+1][1][8][i];
2357	}
2358	L[k+1][2][6][i]=L[k+1][1][7][i];
2359	}
2360	for (k=1; k<=(size(L)+d); k++)
2361	{
2362	L[size(L)][2][5][k]=list();
2363	L[size(L)][2][4][k]=list();
2364	L[size(L)][2][8][k]=list();
2365	L[size(L)][2][3][k]=L[size(L)][2][1][k];
2366	L[size(L)][2][7][k]=L[size(L)][2][6][k];
2367	}
2368	L[size(L)][2][7][size(L)+d+1]=L[size(L)][2][6][size(L)+d+1];
2369	L[size(L)][2][8][size(L)+d+1]=list();
2370	/* building the resolution of the last short exact piece*/
2371	for (i=2; i<=d+size(L); i++)
2372	{
2373	v=0;
2374	if(size(L[size(L)][2][1][i-1])!=0)
2375	{
2376	L[size(L)][2][2][i]=unitmat(nrows(L[size(L)][2][1][i-1]));
2377	}
2378	else
2379	{
2380	L[size(L)][2][2][i-1]=list();
2381	}
2382	}
2383	return(L);
2384	}
2385	/case Syzstring=="Vdres"/
2386	list forVd;
2387	for (k=1; k<=(size(L)+d); k++)//?????
2388	{
2389	/* we compute a V_d-strict resolution for the first short exact piece*/
2390	L[1][1][1][k+1]=list();
2391	L[1][1][2][k+1]=list();
2392	L[1][1][6][k+1]=list();
2393	if (size(L[1][1][3][k])!=0)
2394	{
2395	for (i=1; i<=nrows(L[1][1][3][k]); i++)
2396	{
2397	rem=submat(L[1][1][3][k],i,(1..ncols(L[1][1][3][k])));
2398	n_b[i]=VdDeg(rem,d,L[1][1][7][k]);
2399	}
2400	J_B=transpose(syz(transpose(L[1][1][3][k])));
2401	L[1][1][7][k+1]=n_b;
2402	L[1][1][8][k+1]=n_b;
2403	L[1][1][4][k+1]=unitmat(nrows(L[1][1][3][k]));
2404	if (J_B!=matrix(0,nrows(J_B),ncols(J_B)))
2405	{
2406	J_B=VdStrictGB(J_B,d,n_b);
2407	L[1][1][3][k+1]=J_B;
2408	L[1][1][5][k+1]=J_B;
2409	}
2410	else
2411	{
2412	L[1][1][3][k+1]=list();
2413	L[1][1][5][k+1]=list();
2414	}
2415	n_b=0;
2416	}
2417	else
2418	{
2419	L[1][1][3][k+1]=list();
2420	L[1][1][5][k+1]=list();
2421	L[1][1][7][k+1]=list();
2422	L[1][1][8][k+1]=list();
2423	L[1][1][4][k+1]=list();
2424	}
2425	/* we compute step by step V_d-strict resolutions over
2426	coker(L[i][2][1][1])->coker(L[i][2][3][1])->coker(L[i][2][1][5])
2427	and coker(L[i+1][1][1][1])->coker(L[i+1][1][3][1])->coker(L[i+1][1][1][5])
2428	using the already computed resolutions for coker(L[i][2][1][1])=
2429	coker(L[i][1][3][1]) and coker(L[i+1][1][1][1])=coker(L[i][2][5][1])*/
2430	for (i=1; i<size(L); i++)
2431	{
2432	forVd[1]=L[i][2][1][k];
2433	forVd[2]=L[i][2][2][k];
2434	forVd[3]=L[i][2][3][k];
2435	forVd[4]=L[i][2][4][k];
2436	forVd[5]=L[i][2][5][k];
2437	forVd[6]=L[i][2][6][k];
2438	forVd[7]=L[i][2][7][k];
2439	forVd[8]=L[i][2][8][k];
2440	store=toVdStrict2x3Complex(forVd,d,L[i][1][3][k+1],L[i][1][7][k+1]);
2441	for (j=1; j<=8; j++)
2442	{
2443	L[i][2][j][k+1]=store[j];
2444	}
2445	forVd[1]=L[i+1][1][1][k];
2446	forVd[2]=L[i+1][1][2][k];
2447	forVd[3]=L[i+1][1][3][k];
2448	forVd[4]=L[i+1][1][4][k];
2449	forVd[5]=L[i+1][1][5][k];
2450	forVd[6]=L[i+1][1][6][k];
2451	forVd[7]=L[i+1][1][7][k];
2452	forVd[8]=L[i+1][1][8][k];
2453	store=toVdStrict2x3Complex(forVd,d,L[i][2][5][k+1],L[i][2][8][k+1]);
2454	for (j=1; j<=8; j++)
2455	{
2456	L[i+1][1][j][k+1]=store[j];
2457	}
2458	}
2459	if (size(L[size(L)][1][7][k+1])!=0)
2460	{
2461	L[size(L)][2][4][k+1]=list();
2462	L[size(L)][2][5][k+1]=list();
2463	L[size(L)][2][6][k+1]=L[size(L)][1][7][k+1];
2464	L[size(L)][2][7][k+1]=L[size(L)][1][7][k+1];
2465	L[size(L)][2][8][k+1]=list();
2466	L[size(L)][2][2][k+1]=unitmat(size(L[size(L)][1][7][k+1]));
2467	if (size(L[size(L)][1][3][k+1])!=0)
2468	{
2469	L[size(L)][2][1][k+1]=L[size(L)][1][3][k+1];
2470	L[size(L)][2][3][k+1]=L[size(L)][1][3][k+1];
2471	}
2472	else
2473	{
2474	L[size(L)][2][1][k+1]=list();
2475	L[size(L)][2][3][k+1]=list();
2476	}
2477	}
2478	else
2479	{
2480	for (j=1; j<=8; j++)
2481	{
2482	L[size(L)][2][j][k+1]=list();
2483	}
2484	}
2485	}
2486	k=t;
2487	intvec n_c;
2488	intvec vn_b;
2489	list N_b;
2490	int n;
2491	/computation of the shift vectors/
2492	for (i=1; i<=size(L); i++)
2493	{
2494	for (n=1; n<=2; n++)
2495	{
2496	if (i==1 and n==1)
2497	{
2498	L[i][n][6][k+1]=list();
2499	}
2500	else
2501	{
2502	if (n==1)
2503	{
2504	L[i][1][6][k+1]=L[i-1][2][8][k+1];
2505	}
2506	else
2507	{
2508	L[i][2][6][k+1]=L[i][1][7][k+1];
2509	}
2510	}
2511	N_b[1]=L[i][n][6][k+1];
2512	if (size(L[i][n][5][k])!=0)
2513	{
2514	for (j=1; j<=nrows(L[i][n][5][k]); j++)
2515	{
2516	rem=submat(L[i][n][5][k],j,(1..ncols(L[i][n][5][k])));
2517	n_c[j]=VdDeg(rem,d,L[i][n][8][k]);
2518	}
2519	L[i][n][8][k+1]=n_c;
2520	}
2521	else
2522	{
2523	L[i][n][8][k+1]=list();
2524	}
2525	N_b[2]=L[i][n][8][k+1];
2526	n_c=0;
2527	if (size(N_b[1])!=0)
2528	{
2529	vn_b=N_b[1];
2530	if (size(N_b[2])!=0)
2531	{
2532	vn_b=vn_b,N_b[2];
2533	}
2534	L[i][n][7][k+1]=vn_b;
2535	}
2536	else
2537	{
2538	if (size(N_b[2])!=0)
2539	{
2540	L[i][n][7][k+1]=N_b[2];
2541	}
2542	else
2543	{
2544	L[i][n][7][k+1]=list();
2545	}
2546	}
2547	}
2548	}
2549	return(L);
2550	}
2551
2552	////////////////////////////////////////////////////////////////////////////////////
2553
2554	static proc toVdStrict2x3Complex(list L,int d,list #)
2555	{
2556	/* We build a one-step free resolution over a V_d-strict short exact piece
2557	(Algorithm 3.14 in [W2]).
2558	This procedure is called from the procedure VdStrictDoubleComplexes
2559	if Syzstring=='Vdres'*/
2560	matrix rem;
2561	int i,j,cc;
2562	list J_A=list(list());
2563	list J_B=list(list());
2564	list J_C=list(list());
2565	list g_AB=list(list());
2566	list g_BC=list(list());
2567	list n_a=list(list());
2568	list n_b=list(list());
2569	list n_c=list(list());
2570	intvec n_b1;
2571	matrix fromnf;
2572	intvec i1,i2;
2573	/* compute a one step V_d-strict resolution for L[5]*/
2574	if (size(L[5])!=0)
2575	{
2576	intvec n_c1;
2577	for (i=1; i<=nrows(L[5]); i++)
2578	{
2579	rem=submat(L[5],i,intvec(1..ncols(L[5])));
2580	n_c1[i]=VdDeg(rem,d, L[8]);//new shift vector
2581	}
2582	n_c[1]=n_c1;
2583	J_C[1]=transpose(syz(transpose(L[5])));
2584	if (J_C[1]!=matrix(0,nrows(J_C[1]),ncols(J_C[1])))
2585	{
2586	J_C[1]=VdStrictGB(J_C[1],d,n_c1);
2587	if (size(#[2])!=0)// new shift vector for the resolution of L[1]
2588	{
2589	n_a[1]=#[2];
2590	n_b1=n_a[1],n_c[1];
2591	n_b[1]=n_b1;
2592	matrix zero[nrows(L[1])][nrows(L[5])];
2593	g_AB=concat(unitmat(nrows(L[1])),matrix(0,nrows(L[1]),nrows(L[5])));
2594	if (size(#[1])!=0)
2595	{
2596	J_A=#[1];// one step V_d-strict resolution for L[1]
2597	/* use resolutions of L[1] and L[5] to build a resolution for
2598	L[3]*/
2599	J_B[1]=transpose(matrix(syz(transpose(L[3]))));
2600	matrix P=matrixLift(J_B[1]*prodr(nrows(L[1]),nrows(L[5])),J_C[1]);
2601	matrix Pi[1][ncols(J_B[1])];
2602	matrix Picombined;
2603	for (i=1; i<=nrows(J_C[1]); i++)
2604	{
2605	for (j=1; j<=nrows(J_B[1]);j++)
2606	{
2607	Pi=Pi+P[i,j]*submat(J_B[1],j,intvec(1..ncols(J_B[1])));
2608	}
2609	if (i==1)
2610	{
2611	Picombined=transpose(Pi);
2612	}
2613	else
2614	{
2615	Picombined=concat(Picombined,transpose(Pi));
2616	}
2617	Pi=0;
2618	}
2619	Picombined=transpose(Picombined);
2620	fromnf=VdNormalForm(Picombined,J_A[1],d,n_a[1],n_c[1]);
2621	i1=intvec(1..nrows(Picombined));
2622	i2=intvec((ncols(J_A[1])+1)..ncols(Picombined));
2623	Picombined=concat(fromnf,submat(Picombined,i1,i2));
2624	J_B[1]=transpose(matrix(J_A[1],nrows(J_A[1]),ncols(J_B[1])));
2625	J_B[1]=transpose(concat(J_B[1],transpose(Picombined)));
2626	g_BC=transpose(concat(transpose(zero),unitmat(nrows(L[5]))));
2627	}
2628	else//L[1] is already a resolution
2629	{
2630	//compute a resolution for L[3]
2631	J_B=transpose(matrix(syz(transpose(L[3]))));
2632	matrix P=matrixLift(J_B[1]*prodr(nrows(L[1]),nrows(L[5])),J_C[1]);
2633	matrix Pi[1][ncols(J_B[1])];
2634	matrix Picombined;
2635	for (i=1; i<=nrows(J_C[1]); i++)
2636	{
2637	for (j=1; j<=nrows(J_B[1]);j++)
2638	{
2639	Pi=Pi+P[i,j]*submat(J_B[1],j,intvec(1..ncols(J_B[1])));
2640	}
2641	if (i==1)
2642	{
2643	Picombined=transpose(Pi);
2644	}
2645	else
2646	{
2647	Picombined=concat(Picombined,transpose(Pi));
2648	}
2649	Pi=0;
2650	}
2651	Picombined=transpose(Picombined);
2652	J_B[1]=Picombined;
2653	g_BC=transpose(concat(transpose(zero),unitmat(nrows(L[5]))));
2654	}
2655	}
2656	else
2657	{
2658	n_b=n_c[1];
2659	J_B[1]=J_C[1];
2660	g_BC=unitmat(ncols(J_C[1]));
2661	}
2662	}
2663	else
2664	{
2665	J_C=list(list());// L[5] is already a resolution
2666	if (size(#[2])!=0)
2667	{
2668	matrix zero[nrows(L[1])][nrows(L[5])];
2669	g_BC=transpose(concat(transpose(zero),unitmat(nrows(L[5]))));
2670	n_a[1]=#[2];
2671	n_b1=n_a[1],n_c[1];
2672	n_b[1]=n_b1;
2673	g_AB=concat(unitmat(nrows(L[1])),matrix(0,nrows(L[1]),nrows(L[5])));
2674	if (size(#[1])!=0)
2675	{
2676	J_A=#[1];
2677	/resolution of L[3]/
2678	nr=nrows(J_A[1]);
2679	J_B=concat(J_A[1],matrix(0,nr,nrows(L[3])-nrows(L[1])));
2680	}
2681	}
2682	else
2683	{
2684	n_b=n_c[1];
2685	g_BC=unitmat(ncols(L[5]));
2686	}
2687	}
2688	}
2689	else// L[5]=list();
2690	{
2691	if (size(#[2])!=0)
2692	{
2693	n_a[1]=#[2];
2694	n_b=n_a[1];
2695	g_AB=unitmat(size(n_b[1]));
2696	if (size(#[1])!=0)
2697	{
2698	J_A=#[1];
2699	J_B[1]=J_A[1];// resolution of L[3] equals that of L[1]
2700	}
2701	}
2702	}
2703	list out=(J_A[1],g_AB[1],J_B[1],g_BC[1],J_C[1],n_a[1],n_b[1],n_c[1]);
2704	return (out);
2705	}
2706
2707	////////////////////////////////////////////////////////////////////////////////////
2708
2709	static proc assemblingDoubleComplexes(list L)
2710	{
2711	/* The input is the output of VdStrictDoubleComplexes, we assemble the
2712	resolutions of the L[i][2][3][1] to obtain a V_d-strict free Cartan-Eilenberg
2713	resolution with modules P^i_j (1<=i<=size(L), j>=0) for the seqeunce
2714	coker(L[1][2][3][1])->...->coker(L[size(L)][2][3][1])*/
2715	list out;
2716	int i,j,k,l,oldj,newj,nr,nc;
2717	for (i=1; i<=size(L); i++)
2718	{
2719	out[i]=list(list());
2720	out[i][1][1]=ncols(L[i][2][3][1]);//rank of module P^i_0
2721	if (size(L[i][2][5][1])!=0)
2722	{
2723	/horizontal differential P^i_0->P^(i+1)_0/
2724	nc=ncols(L[i][2][5][1]);
2725	out[i][1][4]=prodr(ncols(L[i][2][3][1])-ncols(L[i][2][5][1]),nc);
2726	}
2727	else
2728	{
2729	/horizontal differential P^i_0->0/
2730	out[i][1][4]=matrix(0,ncols(L[i][2][3][1]),1);
2731	}
2732	oldj=newj;
2733	for (j=1; j<=size(L[i][2][3]);j++)
2734	{
2735	out[i][j][2]=L[i][2][7][j];//shift vector of P^i_{j-1}
2736	if (size(L[i][2][3][j])==0)
2737	{
2738	newj =j;
2739	break;
2740	}
2741	out[i][j+1]=list();
2742	out[i][j+1][1]=nrows(L[i][2][3][j]);//rank of the module P^i_j
2743	out[i][j+1][3]=L[i][2][3][j];//vertical differential P^i_j->P^(i+1)_j
2744	if (size(L[i][2][5][j])!=0)
2745	{
2746	//horizonal differential P^i_j->P^(i-1)_j
2747	nr=nrows(L[i][2][3][j])-nrows(L[i][2][5][j]);
2748	out[i][j+1][4]=(-1)^j*prodr(nr,nrows(L[i][2][5][j]));
2749	}
2750	else
2751	{
2752	/horizontal differential P^i_j->P^(i-1)_j/
2753	out[i][j+1][4]=matrix(0,nrows(L[i][2][3][j]),1);
2754	}
2755	if(j==size(L[i][2][3]))
2756	{
2757	out[i][j+1][2]=L[i][2][7][j+1];//shift vector of P^i_j
2758	newj=j+1;
2759	}
2760	}
2761	if (i>1)
2762	{
2763
2764	for (k=1; k<=Min(list(oldj,newj)); k++)
2765	{
2766	/horizonal differential P^(i-1)_(k-1)->P^i_(k-1)/
2767	nr=nrows(out[i-1][k][4]);
2768	out[i-1][k][4]=matrix(out[i-1][k][4],nr,out[i][k][1]);
2769	}
2770	for (k=newj+1; k<=oldj; k++)
2771	{
2772	/no differential needed/
2773	out[i-1][k]=delete(out[i-1][k],4);
2774	}
2775	}
2776	}
2777	return (out);
2778	}
2779
2780	////////////////////////////////////////////////////////////////////////////////////
2781
2782	static proc totalComplex(list L);
2783	{
2784	/* Input is the output of assemblingDoubleComplexes.
2785	We obtain a complex C^1[m^1]->...->C^(r)[m^r] with differentials d^i and
2786	shift vectors m^i (where C^r is placed in degree size(L)-1).
2787	This complex is dercribed in the list out as follows:
2788	rank(C^i)=out[3i-2]; m_i=out[3i-1] and d^i=out[3i]/
2789	list out;intvec rem1;intvec v; list remsize; int emp;
2790	int i; int j; int c; int d; matrix M; int k; int l;
2791	int n=nvars(basering) div 2;
2792	list K;
2793	for (i=1; i<=n+1; i++)
2794	{
2795	K[i]=list();
2796	}
2797	L=K+L;
2798	for (i=1; i<=size(L); i++)
2799	{
2800	emp=0;
2801	if (size(L[i])!=0)
2802	{
2803	out[3*i-2]=L[i][1][1];
2804	v=L[i][1][1];
2805	rem1=L[i][1][2];
2806	emp=1;
2807	}
2808	else
2809	{
2810	out[3*i-2]=0;
2811	v=0;
2812	}
2813	for (j=i+1; j<=size(L); j++)
2814	{
2815	if (size(L[j])>=j-i+1)
2816	{
2817	out[3i-2]=out[3i-2]+L[j][j-i+1][1];
2818	if (emp==0)
2819	{
2820	rem1=L[j][j-i+1][2];
2821	emp=1;
2822	}
2823	else
2824	{
2825	rem1=rem1,L[j][j-i+1][2];
2826	}
2827	v[size(v)+1]=L[j][j-i+1][1];
2828	}
2829	else
2830	{
2831	v[size(v)+1]=0;
2832	}
2833	}
2834	out[3*i-1]=rem1;
2835	v[size(v)+1]=0;
2836	remsize[i]=v;
2837	}
2838	int o1;
2839	int o2;
2840	for (i=1; i<=size(L)-1; i++)
2841	{
2842	o1=1;
2843	o2=1;
2844	if (size(out[3*i-2])!=0)
2845	{
2846	o1=out[3*i-2];
2847	}
2848	if (size(out[3*i+1])!=0)
2849	{
2850	o2=out[3*i+1];
2851	}
2852	M=matrix(0,o1,o2);
2853	if (size(L[i])!=0)
2854	{
2855	if (size(L[i][1][4])!=0)
2856	{
2857	M=matrix(L[i][1][4],o1,o2);
2858	}
2859	}
2860	c=remsize[i][1];
2861	for (j=i+1; j<=size(L); j++)
2862	{
2863	if (remsize[i][j-i+1]!=0)
2864	{
2865	for (k=c+1; k<=c+remsize[i][j-i+1]; k++)
2866	{
2867	for (l=d+1; l<=d+remsize[i+1][j-i];l++)
2868	{
2869	M[k,l]=L[j][j-i+1][3][(k-c),(l-d)];
2870	}
2871	}
2872	d=d+remsize[i+1][j-i];
2873	if (remsize[i+1][j-i+1]!=0)
2874	{
2875	for (k=c+1; k<=c+remsize[i][j-i+1]; k++)
2876	{
2877	for (l=d+1; l<=d+remsize[i+1][j-i+1];l++)
2878	{
2879	M[k,l]=L[j][j-i+1][4][k-c,l-d];
2880	}
2881	}
2882	c=c+remsize[i][j-i+1];
2883	}
2884	}
2885	else
2886	{
2887	d=d+remsize[i+1][j-i];
2888	}
2889	}
2890	out[3*i]=M;
2891	d=0; c=0;
2892	}
2893	out[3size(L)]=matrix(0,out[3size(L)-2],1);
2894	return (out);
2895
2896	}
2897
2898	////////////////////////////////////////////////////////////////////////////////////
2899	//COMPUTATION OF THE BLOBAL B-FUNCTION
2900	////////////////////////////////////////////////////////////////////////////////////
2901
2902	static proc globalBFun(list L,list #)
2903	{
2904	/*We assume that the basering is the nth Weyl algebra and that L=(L[1],...,L[s]),
2905	where L[i]=(L[i][1],L[i][2]) and L[i][1] is a m_i x n_i-matrix and L[i][2] an
2906	intvec of size n_i.
2907	We compute bounds for the minimal and maximal integer roots of the b-functions
2908	of coker(L[i][1])[L[i][2]], where L[i][2] is the shift vector (cf. Def.
2909	6.1.1 in [R]) by combining Algorithm 6.1.6 in [R] and the method of principal
2910	intersection (cf. Remark 6.1.7 in [R] 2012).
2911	This works ONLY IF ALL B-FUNCTIONS ARE NON-ZERO, but this is the case since this
2912	proc is only called from the procedure deRhamCohomology and the input comes
2913	originally from the procedure toVdstrictFreeComplex*/
2914	if (size(#)==0)//# may contain the Syzstring
2915	{
2916	string Syzstring="Sres";
2917	}
2918	else
2919	{
2920	string Syzstring=#[1];
2921	}
2922	int i,j;
2923	def W=basering;
2924	int n=nvars(W) div 2;
2925	list G0;
2926	ideal I;
2927	for (j=1; j<=size(L); j++)
2928	{
2929	G0[j]=list();
2930	for (i=1; i<=ncols(L[j][1]); i++)
2931	{
2932	G0[j][i]=I;
2933	}
2934	}
2935	list out;
2936	ideal I; poly f;
2937	intvec i1;
2938	for (j=1; j<=size(L); j++)
2939	{
2940	/*if the shift vector L[j][2] is non-zero we have to compute a V_d-strict
2941	Groebner basis of L[j][1] with respect to the zero shift; otherwise L[i][1]
2942	is already a V_d-strict Groebner basis, because it was obtained by the
2943	procedure toVdStrictFreeComplex*/
2944	if (L[j][2]!=intvec(0:size(L[j][2])))
2945	{
2946	if (Syzstring=="Vdres")
2947	{
2948	L[j][1]=VdStrictGB(L[j][1],n);
2949	}
2950	else
2951	{
2952	def HomWeyl=makeHomogenizedWeyl(n);
2953	setring HomWeyl;
2954	list L=fetch(W,L);
2955	L[j][1]=nHomogenize(L[j][1]);
2956	L[j][1]=transpose(matrix(slimgb(transpose(L[j][1]))));
2957	L[j][1]=subst(L[j][1],h,1);
2958	setring W;
2959	L=fetch(HomWeyl,L);
2960	kill HomWeyl;
2961	}
2962	}
2963	for (i=1; i<=ncols(L[j][1]); i++)
2964	{
2965	G0[j][i]=I;
2966	}
2967	for (i=1; i<=nrows(L[j][1]); i++)
2968	{
2969	/*computes the terms of maximal V_d-degree of the biggest non-zero
2970	component of submat(L[j][1],i,(1..ncols(L[j][1])))*/
2971	i1=(1..ncols(L[j][1]));
2972	out=VdDeg(submat(L[j][1],i,i1),n,intvec(0:size(L[j][2])),1);
2973	f=L[j][1][i,out[2]];
2974	G0[j][out[2]]=G0[j][out[2]],f;
2975	G0[j][out[2]]=compress(G0[j][out[2]]);
2976	}
2977	}
2978	list save;
2979	int l;
2980	list weights;
2981	/*bFctIdealModified computes the intersection of G0[j][i] and
2982	x(1)D(1)+...+x(n)D(n) using the method of principal intersection*/
2983	for (j=1; j<=size(G0); j++)
2984	{
2985	for (i=1; i<=size(G0[j]); i++)
2986	{
2987	G0[j][i]=bFctIdealModified(G0[j][i]);
2988	}
2989	for (i=1; i<=size(G0[j]); i++)
2990	{
2991	weights=list();
2992	if (size(G0[j][i])!=0)
2993	{
2994	for (l=i; l<=size(G0[j]); l++)
2995	{
2996	if (size(G0[j][l])!=0)
2997	{
2998	weights[size(weights)+1]=L[j][2][l];
2999	}
3000	}
3001	G0[j][i]=list(G0[j][i][1]+Min(weights),G0[j][i][2]+Max(weights));
3002	}
3003	}
3004	}
3005	list allmin;
3006	list allmax;
3007	for (j=1; j<=size(G0); j++)
3008	{
3009	for (i=1; i<=size(G0[j]); i++)
3010	{
3011	if (size(G0[j][i])!=0)
3012	{
3013	allmin[size(allmin)+1]=G0[j][i][1];
3014	allmax[size(allmax)+1]=G0[j][i][2];
3015	}
3016	}
3017	}
3018	list minmax=list(Min(allmin),Max(allmax));
3019	return(minmax);
3020	}
3021
3022	////////////////////////////////////////////////////////////////////////////////////
3023
3024	static proc exactGlobalBFun(list L,list #)
3025	{
3026	/*We assume that the basering is the nth Weyl algebra and that L=(L[1],...,L[s]),
3027	where L[i]=(L[i][1],L[i][2]) and L[i][1] is a m_i x n_i-matrix and L[i][2] an
3028	intvec of size n_i.
3029	We compute bounds for the minimal and maximal integer roots of the b-functions
3030	of coker(L[i][1])[L[i][2]], where L[i][2] is the shift vector (cf. Def.
3031	6.1.1 in [R]) by combining Algorithm 6.1.6 in [R] and the method of principal
3032	intersection (cf. Remark 6.1.7 in [R] 2012).
3033	This works ONLY IF ALL B-FUNCTIONS ARE NON-ZERO, but this is the case since this
3034	proc is only called from the procedure deRhamCohomology and the input comes
3035	originally from the procedure toVdstrictFreeComplex*/
3036	if (size(#)==0)//# may contain the Syzstring
3037	{
3038	string Syzstring="Sres";
3039	}
3040	else
3041	{
3042	string Syzstring=#[1];
3043	}
3044	int i,j,k;
3045	def W=basering;
3046	int n=nvars(W) div 2;
3047	list G0;
3048	ideal I;
3049	for (j=1; j<=size(L); j++)
3050	{
3051	G0[j]=list();
3052	for (i=1; i<=ncols(L[j][1]); i++)
3053	{
3054	G0[j][i]=I;
3055	}
3056	}
3057	list out;
3058	matrix M;
3059	ideal I; poly f;
3060	intvec i1;
3061	for (j=1; j<=size(L); j++)
3062	{
3063	M=L[j][1];
3064	/*if the shift vector L[j][2] is non-zero we have to compute a V_d-strict
3065	Groebner basis of L[j][1] with respect to the zero shift; otherwise L[i][1]
3066	is already a V_d-strict Groebner basis, because it was obtained by the
3067	procedure toVdStrictFreeComplex*/
3068	for (k=1; k<=ncols(L[j][1]); k++)
3069	{
3070	L[j][1]=permcol(M,1,k);
3071	if (Syzstring=="Vdres")
3072	{
3073	L[j][1]=VdStrictGB(L[j][1],n);
3074	}
3075	else
3076	{
3077	def HomWeyl=makeHomogenizedWeyl(n);
3078	setring HomWeyl;
3079	list L=fetch(W,L);
3080	L[j][1]=nHomogenize(L[j][1]);
3081	L[j][1]=transpose(matrix(slimgb(transpose(L[j][1]))));
3082	L[j][1]=subst(L[j][1],h,1);
3083	setring W;
3084	L=fetch(HomWeyl,L);
3085	kill HomWeyl;
3086	}
3087	for (i=1; i<=nrows(L[j][1]); i++)
3088	{
3089	/*computes the terms of maximal V_d-degree of the biggest non-zero
3090	component of submat(L[j][1],i,(1..ncols(L[j][1])))*/
3091	i1=(1..ncols(L[j][1]));
3092	out=VdDeg(submat(L[j][1],i,i1),n,intvec(0:size(L[j][2])),1);
3093	f=L[j][1][i,out[2]];
3094	if (out[2]==1)
3095	{
3096	G0[j][k]=G0[j][k],f;
3097	G0[j][k]=compress(G0[j][k]);
3098	}
3099	}
3100	}
3101	}
3102	list save;
3103	int l;
3104	list weights;
3105	/*bFctIdealModified computes the intersection of G0[j][i] and
3106	x(1)D(1)+...+x(n)D(n) using the method of principal intersection*/
3107	for (j=1; j<=size(G0); j++)
3108	{
3109	for (i=1; i<=size(G0[j]); i++)
3110	{
3111	G0[j][i]=bFctIdealModified(G0[j][i]);
3112	}
3113	for (i=1; i<=size(G0[j]); i++)
3114	{
3115	if (size(G0[j][i])!=0)
3116	{
3117	G0[j][i]=list(G0[j][i][1]+L[j][2][i],G0[j][i][2]+L[j][2][i]);
3118	}
3119	}
3120	}
3121	list allmin;
3122	list allmax;
3123	for (j=1; j<=size(G0); j++)
3124	{
3125	for (i=1; i<=size(G0[j]); i++)
3126	{
3127	if (size(G0[j][i])!=0)
3128	{
3129	allmin[size(allmin)+1]=G0[j][i][1];
3130	allmax[size(allmax)+1]=G0[j][i][2];
3131	}
3132	}
3133	}
3134	list minmax=list(Min(allmin),Max(allmax));
3135	return(minmax);
3136	}
3137
3138	////////////////////////////////////////////////////////////////////////////////////
3139
3140	static proc bFctIdealModified (ideal I)
3141	{/modified version of the procedure bfunIdeal from bfun.lib/
3142	def B= basering;
3143	int n = nvars(B) div 2;
3144	intvec w=(1:n);
3145	I= initialIdealW(I,-w,w);
3146	poly s; int i;
3147	for (i=1; i<=n; i++)
3148	{
3149	s=s+x(i)*D(i);
3150	}
3151	/pIntersect computes the intersection on s and I/
3152	vector b = pIntersect(s,I);
3153	list RL = ringlist(B); RL = RL[1..4];
3154	RL[2] = list(safeVarName("s"));
3155	RL[3] = list(list("dp",intvec(1)),list("C",intvec(0)));
3156	def @S = ring(RL); setring @S;
3157	vector b = imap(B,b);
3158	poly bs = vec2poly(b);
3159	ring r=0,s,dp;
3160	poly bs=imap(@S,bs);
3161	/find minimal and maximal integer root/
3162	ideal allfac=factorize(bs,1);
3163	list allfacs;
3164	for (i=1; i<=ncols(allfac); i++)
3165	{
3166	allfacs[i]=allfac[i];
3167	}
3168	number testzero;
3169	list zeros;
3170	for (i=1; i<=size(allfacs); i++)
3171	{
3172	if (deg(allfacs[i])==1)
3173	{
3174	testzero=number(subst(allfacs[i],s,0))/leadcoef(allfacs[i]);
3175	if (testzero-int(testzero)==0)
3176	{
3177	zeros[size(zeros)+1]=int(-1)*int(testzero);
3178	}
3179	}
3180	}
3181	if (size(zeros)!=0)
3182	{
3183	list minmax=(Min(zeros),Max(zeros));
3184	}
3185	else
3186	{
3187	list minmax=list();
3188	}
3189	setring B;
3190	return(minmax);
3191	}
3192
3193	////////////////////////////////////////////////////////////////////////////////////
3194
3195	static proc safeVarName (string s)
3196	{/* from the library "bfun.lib"*/
3197	string S = "," + charstr(basering) + "," + varstr(basering) + ",";
3198	s = "," + s + ",";
3199	while (find(S,s) <> 0)
3200	{
3201	s[1] = "@";
3202	s = "," + s;
3203	}
3204	s = s[2..size(s)-1];
3205	return(s)
3206	}
3207
3208	////////////////////////////////////////////////////////////////////////////////////
3209
3210	static proc globalBFunOT(list L,list #)
3211	{
3212	/*this proc is currently not used since globalBFun computes the same output and is
3213	faster, however globalBFun works only for non-zero b-functions!*/
3214	/*We assume that the basering is the nth Weyl algebra and that L=(L[1],...,L[s]),
3215	where L[i]=(L[i][1],L[i][2]) and L[i][1] is a m_i x n_i-matrix and L[i][2] an
3216	intvec of size n_i.
3217	We compute bounds for the minimal and maximal integer roots of the b-functions
3218	of coker(L[i][1])[L[i][2]], where L[i][2] is the shift vector (cf. Def.
3219	6.1.1 in [R]) using Algorithm 6.1.6 in [R].*/
3220	if (size(#)==0)
3221	{
3222	string Syzstring="Sres";
3223	}
3224	else
3225	{
3226	string Syzstring=#[1];
3227	}
3228	int i; int j;
3229	def W=basering;
3230	int n=nvars(W) div 2;
3231	list G0;
3232	ideal I;
3233	intvec i1;
3234	for (j=1; j<=size(L); j++)
3235	{
3236	G0[j]=list();
3237	for (i=1; i<=ncols(L[j][1]); i++)
3238	{
3239	G0[j][i]=I;
3240	}
3241	}
3242	list out;
3243	for (j=1; j<=size(L); j++)
3244	{
3245	if (L[j][2]!=intvec(0:size(L[j][2])))
3246	{
3247	if (Syzstring=="Vdres")
3248	{
3249	L[j][1]=VdStrictGB(L[j][1],n);
3250	}
3251	else
3252	{
3253	def HomWeyl=makeHomogenizedWeyl(n);
3254	setring HomWeyl;
3255	list L=fetch(W,L);
3256	L[j][1]=nHomogenize(L[j][1]);
3257	L[j][1]=transpose(matrix(slimgb(transpose(L[j][1]))));
3258	L[j][1]=subst(L[j][1],h,1);
3259	setring W;
3260	L=fetch(HomWeyl,L);
3261	kill HomWeyl;
3262	}
3263	}
3264	for (i=1; i<=nrows(L[j][1]); i++)
3265	{
3266	i1=(1..ncols(L[j][1]));
3267	out=VdDeg(submat(L[j][1],i,i1),n,intvec(0:size(L[j][2])),1);
3268	G0[j][out[2]][size(G0[j][out[2]])+1]=(out[1]);
3269	}
3270	}
3271	list Data=ringlist(W);
3272	for (i=1; i<=n; i++)
3273	{
3274	Data[2][2*n+i]=Data[2][i];
3275	Data[2][3*n+i]=Data[2][n+i];
3276	Data[2][i]="v("+string(i)+")";
3277	Data[2][n+i]="w("+string(i)+")";
3278	}
3279	Data[3][1][1]="M";
3280	intvec mord=(0:16*n^2);
3281	mord[1..2n]=(1:2n);
3282	mord[6n+1..8n]=(1:2*n);
3283	for (i=0; i<=2*n-2; i++)
3284	{
3285	mord[(3+i)4n-i]=-1;
3286	mord[(2n+2+i)4n-2n-i]=-1;
3287	}
3288	Data[3][1][2]=mord;
3289	matrix Ones=UpOneMatrix(4*n);
3290	Data[5]=Ones;
3291	matrix con[2n][2n];
3292	Data[6]=transpose(concat(con,transpose(concat(con,Data[6]))));
3293	def Wuv=ring(Data);
3294	setring Wuv;
3295	list G0=imap(W,G0); list G3; poly lterm;intvec lexp;
3296	list G1,G2,LL;
3297	intvec e,f;
3298	int kapp,k,l;
3299	poly h;
3300	ideal I;
3301	for (l=1; l<=size(G0); l++)
3302	{
3303	G1[l]=list(); G2[l]=list(); G3[l]=list();
3304	for (i=1; i<=size(G0[l]); i++)
3305	{
3306	for (j=1; j<=ncols(G0[l][i]);j++)
3307	{
3308	G0[l][i][j]=mHom(G0[l][i][j]);
3309	}
3310	for (j=1; j<=nvars(Wuv) div 4; j++)
3311	{
3312	G0[l][i][size(G0[l][i])+1]=1-v(j)*w(j);
3313	}
3314	G1[l][i]=slimgb(G0[l][i]);
3315	G2[l][i]=I;
3316	G3[l][i]=list();
3317	for (j=1; j<=ncols(G1[l][i]); j++)
3318	{
3319	e=leadexp(G1[l][i][j]);
3320	f=e[1..2*n];
3321	if (f==intvec(0:(2*n)))
3322	{
3323	for (k=1; k<=n; k++)
3324	{
3325	kapp=-e[2n+k]+e[3n+k];
3326	if (kapp>0)
3327	{
3328	G1[l][i][j]=(x(k)^kapp)*G1[l][i][j];
3329	}
3330	if (kapp<0)
3331	{
3332	G1[l][i][j]=(D(k)^(-kapp))*G1[l][i][j];
3333	}
3334	}
3335	G2[l][i][size(G2[l][i])+1]=G1[l][i][j];
3336	G3[l][i][size(G3[l][i])+1]=list();
3337	while (G1[l][i][j]!=0)
3338	{
3339	lterm=lead(G1[l][i][j]);
3340	G1[l][i][j]=G1[l][i][j]-lterm;
3341	lexp=leadexp(lterm);
3342	lexp=lexp[2n+1..3n];
3343	LL=list(lexp,leadcoef(lterm));
3344	G3[l][i][size(G3[l][i])][size(G3[l][i][size(G3[l][i])])+1]=LL;
3345	}
3346	}
3347	}
3348	}
3349	}
3350	ring r=0,(s(1..n)),dp;
3351	ideal I;
3352	map G3forr=Wuv,I;
3353	list G3=G3forr(G3);
3354	poly fs,gs;
3355	int a;
3356	list G4;
3357	for (l=1; l<=size(G3); l++)
3358	{
3359	G4[l]=list();
3360	for (i=1; i<=size(G3[l]);i++)
3361	{
3362	G4[l][i]=I;
3363
3364	for (j=1; j<=size(G3[l][i]); j++)
3365	{
3366	fs=0;
3367	for (k=1; k<=size(G3[l][i][j]); k++)
3368	{
3369	gs=1;
3370	for (a=1; a<=n; a++)
3371	{
3372	if (G3[l][i][j][k][1][a]!=0)
3373	{
3374	gs=gs*permuteVar(list(G3[l][i][j][k][1][a]),a);
3375	}
3376	}
3377	gs=gs*G3[l][i][j][k][2];
3378	fs=fs+gs;
3379	}
3380	G4[l][i]=G4[l][i],fs;
3381	}
3382	}
3383	}
3384	if (n==1)
3385	{
3386	ring rnew=0,t,dp;
3387	}
3388	else
3389	{
3390	ring rnew=0,(t,s(2..n)),dp;
3391	}
3392	ideal Iformap;
3393	Iformap[1]=t;
3394	poly forel=1;
3395	for (i=2; i<=n; i++)
3396	{
3397	Iformap[1]=Iformap[1]-s(i);
3398	Iformap[i]=s(i);
3399	forel=forel*s(i);
3400	}
3401	map rtornew=r,Iformap;
3402	list G4=rtornew(G4);
3403	list getintvecs=fetch(W,L);
3404	ideal J;
3405	option(redSB);
3406	for (l=1; l<=size(G4); l++)
3407	{
3408	J=1;
3409	for (i=1; i<=size(G4[l]); i++)
3410	{
3411	G4[l][i]=eliminate(G4[l][i],forel);
3412	J=intersect(J,G4[l][i]);
3413	}
3414	G4[l]=poly(std(J)[1]);
3415	}
3416	list minmax;
3417	list mini=list();
3418	list maxi=list();
3419	list L=fetch(W,L);
3420	for (i=1; i<=size(G4); i++)
3421	{
3422	minmax[i]=minIntRoot0(G4[i],1);
3423	if (size(minmax[i])!=0)
3424	{
3425	mini=insert(mini,minmax[i][1]+Min(L[i][2]));
3426	maxi=insert(maxi,minmax[i][2]+Max(L[i][2]));
3427	}
3428	}
3429	mini=Min(mini);
3430	maxi=Max(maxi);
3431	minmax=list(mini[1],maxi[1]);
3432	option(none);
3433	return(minmax);
3434	}
3435
3436	////////////////////////////////////////////////////////////////////////////////////
3437	//COMPUTATION OF THE COHOMOLOGY
3438	////////////////////////////////////////////////////////////////////////////////////
3439
3440	static proc findCohomology(list L,int le)
3441	{
3442	/computes the cohomology of the complex (D^i,d^i) given by D^i=C^L[2i-1] and
3443	d^i=L[2i]/
3444	def R=basering;
3445	ring r=0,(x),dp;
3446	list L=imap(R,L);
3447	list out;
3448	int i, ker, im;
3449	matrix S;
3450	option(returnSB);
3451	option(redSB);
3452	for (i=2; i<=size(L); i=i+2)
3453	{
3454	if (L[i-1]==0)
3455	{
3456	out[i div 2]=0;
3457	im=0;
3458	}
3459	else
3460	{
3461	S=matrix(syz(transpose(L[i])));
3462	if (S!=matrix(0,nrows(S),ncols(S)))
3463	{
3464	ker=ncols(S);
3465	out[i div 2]=ker-im;
3466	im=L[i-1]-ker;
3467	}
3468	else
3469	{
3470	out[i-1]=0;
3471	im=L[i-1];
3472	}
3473	}
3474	}
3475	option(none);
3476	while (size(out)>le)
3477	{
3478	out=delete(out,1);
3479	}
3480	setring R;
3481	return(out);
3482	}
3483
3484	////////////////////////////////////////////////////////////////////////////////////
3485	//AUXILIARY PROCEDURES
3486	////////////////////////////////////////////////////////////////////////////////////
3487
3488	static proc divdr(matrix m,matrix n)
3489	{
3490	m=transpose(m);
3491	n=transpose(n);
3492	matrix con=concat(m,n);
3493	matrix s=syz(con);
3494	s=submat(s,1..ncols(m),1..ncols(s));
3495	s=transpose(compress(s));
3496	return(s);
3497	}
3498	////////////////////////////////////////////////////////////////////////////////////
3499
3500	static proc matrixLift(matrix M,matrix N)
3501	{
3502	intvec v=option(get);
3503	option(none);
3504	matrix l=transpose(lift(transpose(M),transpose(N)));
3505	option(set,v);
3506	return(l);
3507	}
3508
3509	////////////////////////////////////////////////////////////////////////////////////
3510
3511	static proc VdStrictGB (matrix M,int d,list #)
3512	"USAGE:VdStrictGB(M,d[,v]); M a matrix, d an integer, v an optional intvec
3513	ASSUME:-basering is the nth Weyl algebra D_n @*
3514	-1<=d<=n @*
3515	-v (if given) is the shift vector on the range of M (in particular,
3516	size(v)=ncols(M)); otherwise v is assumed to be the zero shift vector
3517	RETURN:matrix N; the rows of N form a V_d-strict Groebner basis with respect to v
3518	for the module generated by the rows of M
3519	"
3520	{
3521	if (M==matrix(0,nrows(M),ncols(M)))
3522	{
3523	return (matrix(0,1,ncols(M)));
3524	}
3525	intvec op=option(get);
3526	def W =basering;
3527	int ncM=ncols(M);
3528	list Data=ringlist(W);
3529	Data[2]=list("nhv")+Data[2];
3530	Data[3][3]=Data[3][1];
3531	Data[3][1]=list("dp",intvec(1));
3532	matrix re[size(Data[2])][size(Data[2])]=UpOneMatrix(size(Data[2]));
3533	Data[5]=re;
3534	int k,l;
3535	Data[6]=transpose(concat(matrix(0,1,1),transpose(concat(matrix(0,1,1),Data[6]))));
3536	def Whom=ring(Data);// D_n[nhv] with the new commuative variable nhv
3537	setring Whom;
3538	matrix Mnew=imap(W,M);
3539	intvec v;
3540	if (size(#)!=0)
3541	{
3542	v=#[1];
3543	}
3544	if (size(v) < ncM)
3545	{
3546	v=v,0:(ncM-size(v));
3547	}
3548	Mnew=homogenize(Mnew, d, v);//homogenization of M with respect to the new variable
3549	Mnew=transpose(Mnew);
3550	Mnew=slimgb(Mnew);// computes a Groebner basis of the homogenzition of M
3551	Mnew=subst(Mnew,nhv,1);// substitution of 1 gives V_d-strict Groebner basis of M
3552	Mnew=compress(Mnew);
3553	Mnew=transpose(Mnew);
3554	setring W;
3555	M=imap(Whom,Mnew);
3556	option(set,op);
3557	return(M);
3558	}
3559
3560	////////////////////////////////////////////////////////////////////////////////////
3561
3562	static proc VdNormalForm(matrix F,matrix M,int d,intvec v,list #)
3563	"USAGE:VdNormalForm(F,M,d,v[,w]); F and M matrices, d int, v intvec, w an optional
3564	intvec
3565	ASSUME:-basering is the nth Weyl algebra D_n @*
3566	-F a n_1 x n_2-matrix and M a m_1 x m_2-matrix with m_2<=n_2 @*
3567	-d is an integer between 1 and n @*
3568	-v is a shift vector for D_n^(m_2) and hence size(v)=m_2 @*
3569	-w is a shift vector for D_n^(m_1-m_2) and hence size(v)=m_1-m_2 @*
3570	RETURN:a n_1 x n_2-matrix N such that:@*
3571	-If no optional intvec w is given:(N[i,1],..,N[i,m_2]) is a V_d-strict normal
3572	form of (F[i,1],...,F[i,m_2]) with respect to a V_d-strict Groebner basis of
3573	the rows of M and the shift vector v
3574	-If w is given:(N[i,1],..,N[i,m_2]) is chosen such that
3575	Vddeg((N[i,1],...,N[i,m_2])[v])<=Vddeg((F[i,m_2+1],...,F[i,m_1])[v]);
3576	-N[i,j]=F[i,j] for j>m_2
3577	"
3578	{
3579	int SBcom;
3580	def W =basering;
3581	int c=ncols(M);
3582	matrix keepF=F;
3583	if (size(#)!=0)
3584	{
3585	intvec w=#[1];
3586	}
3587	F=submat(F,intvec(1..nrows(F)),intvec(1..c));
3588	list Data=ringlist(W);
3589	Data[2]=list("nhv")+Data[2];
3590	Data[3][3]=Data[3][1];
3591	Data[3][1]=list("dp",intvec(1));
3592	matrix re[size(Data[2])][size(Data[2])]=UpOneMatrix(size(Data[2]));
3593	Data[5]=re;
3594	int k,l,nr,nc;
3595	matrix rep[size(Data[2])][size(Data[2])];
3596	for (l=size(Data[2])-1;l>=1; l--)
3597	{
3598	for (k=l-1; k>=1;k--)
3599	{
3600	rep[k+1,l+1]=Data[6][k,l];
3601	}
3602	}
3603	Data[6]=rep;
3604	def Whom=ring(Data);//new ring D_n[nvh] this new commuative variable nhv
3605	setring Whom;
3606	matrix Mnew=imap(W,M);
3607	list forMnew=homogenize(Mnew,d,v,1);//commputes homogenization of M;
3608	Mnew=forMnew[1];
3609	int rightexp=forMnew[2];
3610	matrix Fnew=imap(W,F);
3611	matrix keepF=imap(W,keepF);
3612	matrix Fb;
3613	int cc;
3614	intvec i1,i2;
3615	matrix zeromat,subm1,subm2,zeromat2;
3616	for (l=1; l<=nrows(Fnew); l++)
3617	{
3618	if (size(#)!=0)
3619	{
3620	subm2=submat(keepF,l,((ncols(Fnew)+1)..ncols(keepF)));
3621	zeromat2=matrix(0,1,ncols(subm2));
3622	if (submat(keepF,l,((ncols(Fnew)+1)..ncols(keepF)))==zeromat2)
3623	{
3624	for (cc=1; cc<=ncols(Fnew); c++)
3625	{
3626	Fnew[l,cc]=0;
3627	}
3628	}
3629	i1=intvec(1..ncols(Fnew));
3630	subm1=submat(Fnew,l,i1);
3631	subm2=submat(keepF,l,(ncols(Fnew)+1)..ncols(keepF));
3632	zeromat=matrix(0,1,ncols(Fnew));
3633	if (VdDegnhv(subm1,d,v)>VdDegnhv(subm2,d,w)
3634	and submat(Fnew,l,intvec(1..ncols(Fnew)))!=zeromat)
3635	{
3636	//print("Reduzierung des V_d-Grades nötig");
3637	/*We need to reduce the V_d-degree. First we homogenize the
3638	lth row of Fnew*/
3639	Fb=homogenize(subm1,d,v)*(nhv^rightexp);
3640	if (SBcom==0)
3641	{
3642	/computes a V_d-strict standard basis/
3643	Mnew=slimgb(transpose(Mnew));//
3644	SBcom=1;
3645	}
3646	/computes a V_d-strict normal form for FB/
3647	Fb=transpose(reduce(transpose(Fb),groebner(Mnew)));
3648	if (VdDegnhv(Fb,d,v)> VdDegnhv(subm2,d,w)
3649	and Fb!=matrix(0,nrows(Fb),ncols(Fb)))//should not happen
3650	{
3651	//print("Reduzierung fehlgeschlagen!!!!!!!!!!!!!!!!");
3652	}
3653	}
3654	else
3655	{
3656	/*condition on V_ddeg already satisfied -> no normal form
3657	computation is needed*/
3658	Fb=submat(Fnew,l,intvec(1..ncols(Fnew)));
3659	}
3660	}
3661	else
3662	{
3663	Fb=homogenize(submat(Fnew,l,intvec(1..ncols(Fnew))),d,v);
3664	if (SBcom==0)
3665	{
3666	Mnew=slimgb(transpose(Mnew));// computes a V_d-strict Groebner basis
3667	SBcom=1;
3668	}
3669	Fb=transpose(reduce(transpose(Fb),groebner(Mnew)));//normal form
3670	}
3671	for (k=1; k<=ncols(Fnew);k++)
3672	{
3673	Fnew[l,k]=Fb[1,k];
3674	}
3675	}
3676	Fnew=subst(Fnew,nhv,1);//obtain normal form in D_n
3677	setring W;
3678	F=imap(Whom,Fnew);
3679	return(F);
3680	}
3681
3682	////////////////////////////////////////////////////////////////////////////////////
3683
3684	static proc homogenize (matrix M,int d,intvec v,list #)
3685	{
3686	/* we compute the F[v]-homogenization of each row of M (cf. Def. 3.4 in [OT])*/
3687	if (M==matrix(0,nrows(M),ncols(M)))
3688	{
3689	return(M);
3690	}
3691	int i,l,s, kmin, nhvexp;
3692	poly f;
3693	intvec vnm;
3694	list findmin,maxnhv,rempoly,remk,rem1,rem2;
3695	int n=(nvars(basering)-1) div 2;
3696	for (int k=1; k<=nrows(M); k++)
3697	{
3698	for (l=1; l<=ncols (M); l++)
3699	{
3700	f=M[k,l];
3701	s=size(f);
3702	for (i=1; i<=s; i++)
3703	{
3704	vnm=leadexp(f);
3705	vnm=vnm[n+2..n+d+1]-vnm[2..d+1];
3706	kmin=sum(vnm)+v[l];
3707	rem1[size(rem1)+1]=lead(f);
3708	rem2[size(rem2)+1]=kmin;
3709	findmin=insert(findmin,kmin);
3710	f=f-lead(f);
3711	}
3712	rempoly[l]=rem1;
3713	remk[l]=rem2;
3714	rem1=list();
3715	rem2=list();
3716	}
3717	if (size(findmin)!=0)
3718	{
3719	kmin=Min(findmin);
3720	}
3721	for (l=1; l<=ncols(M); l++)
3722	{
3723	if (M[k,l]!=0)
3724	{
3725	M[k,l]=0;
3726	for (i=1; i<=size(rempoly[l]);i++)
3727	{
3728	nhvexp=remk[l][i]-kmin;
3729	M[k,l]=M[k,l]+nhv^(nhvexp)*rempoly[l][i];
3730	maxnhv[size(maxnhv)+1]=nhvexp;
3731	}
3732	}
3733	}
3734	rempoly=list();
3735	remk=list();
3736	findmin=list();
3737	}
3738	maxnhv=Max(maxnhv);
3739	nhvexp=maxnhv[1];
3740	if (size(#)!=0)
3741	{
3742	return(list(M,nhvexp));//only needed for normal form computations
3743	}
3744	return(M);
3745	}
3746
3747	////////////////////////////////////////////////////////////////////////////////////
3748
3749	static proc soldr (matrix M,matrix N)
3750	{
3751	/* We compute a ncols(M) x nrows(M)-matrix C such that
3752	C[i,1]M_1+...+C[i,nrows(M)]M_(nrows(M))= e_i mod im(N),
3753	where e_i is the ith basis element on the range of M, M_j denotes the jth row
3754	of M and im(N) is generated by the rows of N */
3755	int n=nrows(M);
3756	int q=ncols(M);
3757	matrix S=concat(transpose(M),transpose(N));
3758	def W=basering;
3759	list Data=ringlist(W);
3760	list Save=Data[3];
3761	Data[3]=list(list("c",0),list("dp",intvec(1..nvars(W))));
3762	def Wmod=ring(Data);
3763	setring Wmod;
3764	matrix Smod=imap(W,S);
3765	matrix E[q][1];
3766	matrix Smod2,Smodnew;
3767	option(returnSB);
3768	int i,j;
3769	for (i=1;i<=q;i++)
3770	{
3771	E[i,1]=1;
3772	Smod2=concat(E,Smod);
3773	Smod2=syz(Smod2);
3774	E[i,1]=0;
3775	for (j=1;j<=ncols(Smod2);j++)
3776	{
3777	if (Smod2[1,j]==1)
3778	{
3779	Smodnew=concat(Smodnew,(-1)*(submat(Smod2,intvec(2..n+1),j)));
3780	break;
3781	}
3782	}
3783	}
3784	Smodnew=transpose(submat(Smodnew,intvec(1..n),intvec(2..q+1)));
3785	setring W;
3786	matrix Snew=imap(Wmod,Smodnew);
3787	option(none);
3788	return (Snew);
3789	}
3790
3791	////////////////////////////////////////////////////////////////////////////////////
3792
3793	static proc prodr (int k,int l)
3794	{
3795	if (k==0)
3796	{
3797	matrix P=unitmat(l);
3798	return (P);
3799	}
3800	matrix O[l][k];
3801	matrix P=transpose(concat(O,unitmat(l)));
3802	return (P);
3803	}
3804
3805	////////////////////////////////////////////////////////////////////////////////////
3806
3807	static proc VdDeg(matrix M,int d,intvec v,list #)
3808	{
3809	/* We assume that the basering it the nth Weyl algebra and that M is a 1 x r-
3810	matrix.
3811	We compute the V_d-deg of M with respect to the shift vector v,
3812	i.e V_ddeg(M)=max (V_ddeg(M_i)+v[i]), where k=V_ddeg(M_i) if k is the minimal
3813	integer, such that M_i can be expressed as a sum of operators
3814	x(1)^(a_1)...x(n)^(a_n)D(1)^(b_1)...*D(n)^(b_n) with
3815	a_1+..+a_d+k>=b_1+..+b_d*/
3816	int i, j, etoint;
3817	int n=nvars(basering) div 2;
3818	intvec e;
3819	list findmax;
3820	int c=ncols(M);
3821	poly l;
3822	list positionpoly,positionVd;
3823	for (i=1; i<=c; i++)
3824	{
3825	positionpoly[i]=list();
3826	positionVd[i]=list();
3827	while (M[1,i]!=0)
3828	{
3829	l=lead(M[1,i]);
3830	positionpoly[i][size(positionpoly[i])+1]=l;
3831	e=leadexp(l);
3832	e=-e[1..d]+e[n+1..n+d];
3833	e=sum(e)+v[i];
3834	etoint=e[1];
3835	positionVd[i][size(positionVd[i])+1]=etoint;
3836	findmax[size(findmax)+1]=etoint;
3837	M[1,i]=M[1,i]-l;
3838	}
3839	}
3840	if (size(findmax)!=0)
3841	{
3842	int maxVd=Max(findmax);
3843	if (size(#)==0)
3844	{
3845	return (maxVd);
3846	}
3847	}
3848	else // M is 0-modul
3849	{
3850	return(int(0));
3851	}
3852	l=0;
3853	for (i=c; i>=1; i--)
3854	{
3855	for (j=1; j<=size(positionVd[i]); j++)
3856	{
3857	if (positionVd[i][j]==maxVd)
3858	{
3859	l=l+positionpoly[i][j];
3860	}
3861	}
3862	if (l!=0)
3863	{
3864	/*returns the largest component that has maximal V_d-degree
3865	and its terms of maximal V_d-deg (needed for globalBFun)*/
3866	return (list(l,i));
3867	}
3868	}
3869	}
3870
3871	////////////////////////////////////////////////////////////////////////////////////
3872
3873	static proc VdDegnhv(matrix M,int d,intvec v,list #)
3874	{
3875	/* As the procedure VdDeg, but the basering is the nth Weyl algebra
3876	with a commutative variable nhv*/
3877	int i,j,etoint;
3878	int n=nvars(basering) div 2;
3879	intvec e;
3880	int etoint;
3881	list findmax;
3882	int c=ncols(M);
3883	poly l;
3884	list positionpoly;
3885	list positionVd;
3886	for (i=1; i<=c; i++)
3887	{
3888	positionpoly[i]=list();
3889	positionVd[i]=list();
3890	while (M[1,i]!=0)
3891	{
3892	l=lead(M[1,i]);
3893	positionpoly[i][size(positionpoly[i])+1]=l;
3894	e=leadexp(l);
3895	e=-e[2..d+1]+e[n+2..n+d+1];
3896	e=sum(e)+v[i];
3897	etoint=e[1];
3898	positionVd[i][size(positionVd[i])+1]=etoint;
3899	findmax[size(findmax)+1]=etoint;
3900	M[1,i]=M[1,i]-l;
3901	}
3902	}
3903	if (size(findmax)!=0)
3904	{
3905	int maxVd=Max(findmax);
3906	if (size(#)==0)
3907	{
3908	return (maxVd);
3909	}
3910	}
3911	else // M is 0-modul
3912	{
3913	return(int(0));
3914	}
3915	}
3916
3917	////////////////////////////////////////////////////////////////////////////////////
3918
3919	static proc deletecol(matrix M,int l)
3920	{
3921	int s=ncols(M);
3922	if (l==1)
3923	{
3924	M=submat(M,(1..nrows(M)),(2..ncols(M)));
3925	return(M);
3926	}
3927	if (l==s)
3928	{
3929	M=submat(M,(1..nrows(M)),(1..(ncols(M)-1)));
3930	return(M);
3931	}
3932	intvec v=(1..(l-1)),((l+1)..s);
3933	M=submat(M,(1..nrows(M)),v);
3934	return(M);
3935	}
3936
3937	////////////////////////////////////////////////////////////////////////////////////
3938
3939	static proc mHom(poly f)
3940	{/for globalBFunOT/
3941	poly g;
3942	poly l;
3943	poly add;
3944	intvec e;
3945	list minint;
3946	list remf;
3947	int i;
3948	int j;
3949	int n=nvars(basering) div 4;
3950	if (f==0)
3951	{
3952	return(f);
3953	}
3954	while (f!=0)
3955	{
3956	l=lead(f);
3957	e=leadexp(l);
3958	remf[size(remf)+1]=list();
3959	remf[size(remf)][1]=l;
3960	for (i=1; i<=n; i++)
3961	{
3962	remf[size(remf)][i+1]=-e[2n+i]+e[3n+i];
3963	if (size(minint)<i)
3964	{
3965	minint[i]=list();
3966	}
3967	minint[i][size(minint[i])+1]=-e[2n+i]+e[3n+i];
3968	}
3969	f=f-l;
3970	}
3971	for (i=1; i<=n; i++)
3972	{
3973	minint[i]=Min(minint[i]);
3974	}
3975	for (i=1; i<=size(remf); i++)
3976	{
3977	add=remf[i][1];
3978	for (j=1; j<=n; j++)
3979	{
3980	add=v(j)^(remf[i][j+1]-minint[j])*add;
3981	}
3982	g=g+add;
3983	}
3984	return (g);
3985	}
3986
3987	////////////////////////////////////////////////////////////////////////////////////
3988
3989	static proc permuteVar(list L,int n)
3990	{/for globalBFunOT/
3991	if (typeof(L[1])=="intvec")
3992	{
3993	intvec v=L[1];
3994	}
3995	else
3996	{
3997	intvec v=(1:L[1]),(0:L[1]);
3998	}
3999	int i;int k; int indi=0;
4000	int j;
4001	int s=size(v);
4002	poly e;
4003	intvec fore;
4004	for (i=2; i<=size(v); i=i+2)
4005	{
4006
4007	if (v[i]!=0)
4008	{
4009	j=i+1;
4010	while (v[j]!=0)
4011	{
4012	j=j+1;
4013	}
4014	v[i]=0;
4015	v[j]=1;
4016	fore=0;
4017	indi=0;
4018	for (k=1; k<=size(v); k++)
4019	{
4020	if (k!=i and k!=j)
4021	{
4022	if (indi==0)
4023	{
4024	indi=1;
4025	fore[1]=v[k];
4026	}
4027	else
4028	{
4029	fore[size(fore)+1]=v[k];
4030	}
4031	}
4032	}
4033	e=e-(j-i)*permutevar(list(fore),n);
4034	}
4035	}
4036	e=e+s(n)^(size(v) div 2);
4037	return (e);
4038	}
4039
4040	////////////////////////////////////////////////////////////////////////////////////
4041
4042	static proc makeHomogenizedWeyl(int n,list #)
4043	{
4044	/modified version of the procedure makeWeyl() from the library nctools.lib/
4045	/*Creates the nth homogenized Weyl algebra with variables x(1),..,x(n),D(1),..,
4046	D(n) and homogenization variable h, i.e. it holds x(i)D(i)=D(i)x(1)+h^2.
4047	If # contains on intvec v, we assign weight v[i] to the ith module component.*/
4048	if (n<1)
4049	{
4050	print("Incorrect input");
4051	return();
4052	}
4053	if (n ==1)
4054	{
4055	ring @rr = 0,(x(1),D(1),h),dp;
4056	}
4057	else
4058	{
4059	ring @rr = 0,(x(1..n),D(1..n),h),dp;
4060	}
4061	setring @rr;
4062	if (size(#)==0)
4063	{
4064	def @rrr = homogenizedWeyl();
4065	}
4066	else
4067	{
4068	def @rrr=homogenizedWeyl(#);
4069	}
4070	return(@rrr);
4071	}
4072
4073	////////////////////////////////////////////////////////////////////////////////////
4074
4075	static proc homogenizedWeyl (list #)
4076	{
4077	/modified version of the procedure Weyl() from the library nctools.lib/
4078	/*Creates a homogenized Weyl algebra structure on the basering. We assume
4079	n=nvars(basering) is odd. The first (n-1)/2 variables will be treated as the
4080	x(i), the next (n-1)/2 as the corresponding differentials D(i) and the last as
4081	the homogenization variable h, i.e. it holds x(i)D(i)=D(i)x(1)+h^2.
4082	If # contains on intvec v, we assign weight v[i] to the ith module component.*/
4083	string rname=nameof(basering);
4084	if ( rname == "basering") // i.e. no ring has been set yet
4085	{
4086	"You have to call the procedure from the ring";
4087	return();
4088	}
4089	int nv = nvars(basering);
4090	int N = (nv-1) div 2;
4091	if (((nv-1) % 2) != 0)
4092	{
4093	"Cannot create homogenized Weyl structure for an even number of generators";
4094	return();
4095	}
4096	matrix @D[nv][nv];
4097	int i;
4098	for ( i=1; i<=N; i++ )
4099	{
4100	@D[i,N+i]=h^2;
4101	}
4102	def @R = nc_algebra(1,@D);
4103	setring @R;
4104	list RL=ringlist(@R);
4105	intvec v;
4106	/we need this ordering for Groebner basis computations/
4107	for (i=1; i<=N; i++)
4108	{
4109	v[i]=-1;
4110	v[N+i]=1;
4111	}
4112	v[nv]=0;
4113	/* we assign weights to module components*/
4114	if (size(#)!=0)
4115	{
4116	if (typeof(#[1])=="intvec")
4117	{
4118	intvec m=#[1];
4119	for (i=1; i<=size(m); i++)
4120	{
4121	v[size(v)+1]=m[i];//assigns weight m[i] to the ith module component
4122	}
4123	RL[3]=insert(RL[3],list("am",v));
4124	}
4125	else
4126	{
4127	RL[3]=insert(RL[3],list("a",v));
4128	}
4129	}
4130	else
4131	{
4132	RL[3]=insert(RL[3],list("a",v));
4133	}
4134	intvec w=(1:nv);
4135	if (size(#)>=2)
4136	{
4137	if (typeof(#[2])=="intvec")
4138	{
4139	intvec n=#[2];
4140	for (i=1; i<=size(n); i++)
4141	{
4142	w[size(w)+1]=n[i];
4143	}
4144	RL[3]=insert(RL[3],list("am",w));
4145	}
4146	else
4147	{
4148	RL[3]=insert(RL[3],list("a",w));
4149	}
4150	}
4151	else
4152	{
4153	RL[3]=insert(RL[3],list("a",w));
4154	}
4155	/this ordering is needed for globalBFun and globalBFunOT/
4156	list saveord=RL[3][3];
4157	RL[3][3]=RL[3][4];
4158	RL[3][4]=saveord;
4159	intvec notforh=(1:(size(RL[3][4][2])-1));
4160	RL[3][4][2]=notforh;
4161	RL[3][5]=list("dp",1);
4162	def @@R=ring(RL);
4163	return(@@R);
4164	}
4165
4166	////////////////////////////////////////////////////////////////////////////////////
4167
4168	static proc nHomogenize (matrix M,list #)
4169	{
4170	/* # may contain an intvec v, if no intvec is given, we assume that v=(0:ncols(M))
4171	We compute the h[v]-homogenization of the rows of M as in Definition 9.2 [OT]*/
4172	int l; poly f; int s; int i; intvec vnm;int kmin; list findmax;
4173	int n=(nvars(basering)-1) div 2;
4174	list rempoly;
4175	list remk;
4176	list rem1;
4177	list rem2;
4178	list maxhexp;
4179	int hexp;
4180	intvec v=(0:ncols(M));
4181	if (size(#)!=0)
4182	{
4183	if (typeof(#[1])=="intvec")
4184	{
4185	v=#[1];
4186	}
4187	}
4188	if (size(v)<ncols(M))
4189	{
4190	for (i=size(v)+1; i<=ncols(M); i++)
4191	{
4192	v[i]=0;
4193	}
4194	}
4195	for (int k=1; k<=nrows(M); k++)
4196	{
4197	for (l=1; l<=ncols (M); l++)
4198	{
4199	f=M[k,l];
4200	s=size(f);
4201	for (i=1; i<=s; i++)
4202	{
4203	vnm=leadexp(f);
4204	kmin=sum(vnm)+v[l];
4205	rem1[size(rem1)+1]=lead(f);
4206	rem2[size(rem2)+1]=kmin;
4207	findmax=insert(findmax,kmin);
4208	f=f-lead(f);
4209	}
4210	rempoly[l]=rem1;
4211	remk[l]=rem2;
4212	rem1=list();
4213	rem2=list();
4214	}
4215	if (size(findmax)!=0)
4216	{
4217	kmin=Max(findmax);
4218	}
4219	else
4220	{
4221	kmin=0;
4222	}
4223	for (l=1; l<=ncols(M); l++)
4224	{
4225	if (M[k,l]!=0)
4226	{
4227	M[k,l]=0;
4228	for (i=1; i<=size(rempoly[l]);i++)
4229	{
4230	hexp=kmin-remk[l][i];
4231	maxhexp[size(maxhexp)+1]=hexp;
4232	M[k,l]=M[k,l]+h^hexp*rempoly[l][i];
4233	}
4234	}
4235	}
4236	rempoly=list();
4237	remk=list();
4238	findmax=list();
4239	}
4240	if (size(maxhexp)!=0)
4241	{
4242	maxhexp=Max(maxhexp);
4243	hexp=maxhexp[1];
4244	}
4245	else
4246	{
4247	hexp=0;
4248	}
4249	if (size(#)>1)
4250	{
4251	list forreturn=M,hexp;
4252
4253	return(forreturn);
4254	}
4255	return(M);
4256	}
4257
4258	////////////////////////////////////////////////////////////////////////////////////
4259
4260	static proc max(int i,int j)
4261	{
4262	if(i>j){return(i);}
4263	return(j);
4264	}
4265
4266	////////////////////////////////////////////////////////////////////////////////////
4267
4268	static proc nDeg (matrix M,intvec m)
4269	{/*we compute an intvec n such that n[i]=max(deg(M[i,j])+m[j]\|M[i,j]!=0) (where deg
4270	stands for the total degree) if (M[i,j]!=0 for some j) and n[i]=0 else*/
4271	int i; int j;
4272	intvec n;
4273	list L;
4274	for (i=1; i<=nrows(M); i++)
4275	{
4276	L=list();
4277	for (j=1; j<=ncols(M); j++)
4278	{
4279	if (M[i,j]!=0)
4280	{
4281	L=insert(L,deg(M[i,j])+m[j]);
4282	}
4283	}
4284	if (size(L)==0)
4285	{
4286	n[i]=0;
4287	}
4288	else
4289	{
4290	n[i]=Max(L);
4291	}
4292	}
4293	return(n);
4294	}
4295
4296	////////////////////////////////////////////////////////////////////////////////////
4297
4298	static proc minIntRoot0(list L,list #)
4299	"USAGE:minIntRoot0(L [,M]); L list, M optinonal list
4300	ASSUME:L a list of univariate polynomials with rational coefficients @*
4301	the variable of the polynomial is s if size(#)==0 (needed for proc
4302	MVComplex) and t else (needed for globalBFun)
4303	RETURN:-if size(#)==0: int i, where i is an integer root of one of the polynomials
4304	and it is minimal with respect to that property@*
4305	-if size(#)!=0: list L=(i,j), where i is as above and j is an integer root
4306	of one of the polynomials and is maximal with respect to that property (if
4307	an integer root exists) or L=list() else
4308	"
4309	{
4310	def B=basering;
4311	if (size(#)==0)
4312	{
4313	ring rnew=0,s,dp;
4314	}
4315	else
4316	{
4317	ring rnew=0,t,dp;
4318	}
4319	list L=imap(B,L);
4320
4321	int i;
4322	int j;
4323	number isint;
4324	list possmin;
4325	ideal allfac;
4326	list allfacs;
4327	for (i=1; i<=size(L); i++)
4328	{
4329	allfac=factorize(L[i],1);
4330	for (j=1; j<=ncols(allfac); j++)
4331	{
4332	allfacs[j]=allfac[j];
4333	}
4334	for (j=1; j<=size(allfacs); j++)
4335	{
4336	if (deg(allfacs[j])==1)
4337	{
4338	isint=number(subst(allfacs[j],var(1),0)/leadcoef(allfacs[j]));
4339	if (isint-int(isint)==0)
4340	{
4341	possmin[size(possmin)+1]=int(isint);
4342	}
4343	}
4344	}
4345	allfacs=list();
4346	}
4347	int zerolist;
4348	if (size(possmin)!=0)
4349	{
4350	int miniroot=(-1)*Max(possmin);
4351	int maxiroot=(-1)*Min(possmin);
4352	}
4353	else
4354	{
4355	zerolist=1;
4356	}
4357	setring B;
4358	if (size(#)==0)
4359	{
4360	return(miniroot);
4361	}
4362	else
4363	{
4364	if (zerolist==0)
4365	{
4366	return(list(miniroot,maxiroot));
4367	}
4368	else
4369	{
4370	return(list());
4371	}
4372	}
4373	}
4374
4375
4376	////////////////////////////////////////////////////////////////////////////////////
4377	////////////////////////////////////////////////////////////////////////////////////
4378	////////////////////////////////////////////////////////////////////////////////////
4379	/*
4380	////////////////////////////////////////////////////////////////////////////////////
4381	FURTHER EXAMPLES FOR TESTING THE PROCEDURES
4382	////////////////////////////////////////////////////////////////////////////////////
4383	LIB "derham.lib";
4384
4385	//----------------------------------------
4386	//EXAMPLE 1
4387	//----------------------------------------
4388	ring r=0,(x,y),dp;
4389	poly f=y2-x3-2x+3;
4390	list L=deRhamCohomology(f);
4391	L;
4392	kill r;
4393
4394	//----------------------------------------
4395	//EXAMPLE 2
4396	//----------------------------------------
4397	ring r=0,(x,y),dp;
4398	poly f=y2-x3-x;
4399	list L=deRhamCohomology(f);
4400	L;
4401	kill r;
4402
4403	//----------------------------------------
4404	//EXAMPLE 3
4405	//----------------------------------------
4406	ring r=0,(x,y),dp;
4407	list C=list(x2-1,(x+1)y,y(y2+2x+1));
4408	list L=deRhamCohomology(C);
4409	L;
4410	kill r;
4411
4412	//----------------------------------------
4413	//EXAMPLE 4
4414	//----------------------------------------
4415	ring r=0,(x,y,z),dp;
4416	list C=list(x(x-1),y,z(z-1),z*(x-1));
4417	list L=deRhamCohomology(C);
4418	L;
4419	kill r;
4420
4421	//----------------------------------------
4422	//EXAMPLE 5
4423	//----------------------------------------
4424	ring r=0,(x,y,z),dp;
4425	list C=list(xy,yz);
4426	list L=deRhamCohomology(C,"Vdres");
4427	L;
4428	kill r;
4429
4430	//----------------------------------------
4431	//EXAMPLE 6
4432	//----------------------------------------
4433	ring r=0,(x,y,z,u),dp;
4434	list C=list(x,y,z,u);
4435	list L=deRhamCohomology(C);
4436	L;
4437	kill r;
4438
4439	//----------------------------------------
4440	//EXAMPLE 7
4441	//----------------------------------------
4442	ring r=0,(x,y,z),dp;
4443	poly f=x3+y3+z3;
4444	list L=deRhamCohomology(f);
4445	L;
4446	kill r;
4447
4448	//----------------------------------------
4449	//EXAMPLE 8
4450	//----------------------------------------
4451	ring r=0,(x,y,z),dp;
4452	poly f=x2+y2+z2;
4453	list L=deRhamCohomology(f,"Vdres");
4454	L;
4455	kill r;
4456
4457	//----------------------------------------
4458	//EXAMPLE 9
4459	//----------------------------------------
4460	ring r=0,(x,y,z,u),dp;
4461	list C=list(x2+y2+z2,u);
4462	list L=deRhamCohomology(C);
4463	L;
4464	kill r;
4465
4466
4467	//----------------------------------------
4468	//EXAMPLE 10
4469	//----------------------------------------
4470	ring r=0,(x,y,z),dp;
4471	list C=list((x*(y-1),y2-1));
4472	list L=deRhamCohomology(C);
4473	L;
4474	kill r;
4475
4476
4477	*/

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

Download in other formats:

Original Format