Context Navigation

invar.lib @ 4b56a8

spielwiese

Last change on this file since 4b56a8 was 4b56a8, checked in by Hans Schoenemann <hannes@…>, 4 years ago
fix: avoid warning in invar.lib
Property mode set to `100644`
File size: 17.8 KB

Line
1	// Harm Derksen, hderksen@math.unibas.ch
2	////////////////////////////////////////////////////////////////////////////////
3	version="version invar.lib 4.1.2.0 Feb_2019 "; // $Id$
4	category="Invariant theory";
5	info="
6	LIBRARY: invar.lib Procedures to compute invariants ring of SL(n) and torus groups
7	AUTHOR: Harm Derksen, hderksen@math.unibas.ch
8
9	PROCEDURES:
10	SL(n) sets the current group to SL_n
11	torus(n) sets the current group to an n-dimensional torus
12	torusrep(list m) representation of a torus given by the weights m[1],m[2],...
13	finiterep(<list>) representation of a by a list of matrices
14	sympower(m,d) computes the d-th symmetric power of a representation m
15	invar(m) computes the invariant ring of the representation m.
16	SLreynolds(f) applies the Reynolds operator to f
17	torusreynolds(f) applies the Reynolds operator to f if the group is a torus or a finite group.
18	cyclotomic(n,z) the n-th cyclotomic polynomial in the variable z.
19	";
20
21	////////////////////////////////////////////////////////////////////////////////
22	// SL(n) sets the current group to SL_n.
23	////////////////////////////////////////////////////////////////////////////////
24	proc SL(int n)
25	"USAGE: SL(<int>)
26	RETURNS: SL(n) sets the current group to SL_n. The following global variables
27	will be set:
28	group of type <ring>
29	groupideal of type <ideal>
30	SLrep of type <matrix>
31	reynolds of type <proc>
32	The quotient of of `group` and `groupideal` is the coordinate
33	ring of SL_n. The matrix `SLrep` will be set to the standard
34	representation of SL_n. The basering will be set to 'group'.
35	EXAMPLE: example SL; shows an example"
36	{ ring group=0,(g(1..n^2)),dp;
37	export group;
38	matrix SLrep[n][n]; // SLrep is a generic n x n
39	int i; // matrix
40	int j;
41	for(i=1;i<=n;i=i+1)
42	{
43	for(j=1;j<=n;j=j+1)
44	{SLrep[i,j]=g(j+n*(i-1));}
45	}
46	ideal groupideal=det(SLrep)-1; // SL(n) are those matrices m
47	groupideal=std(groupideal); // with determinant equal to 1.
48	export groupideal;
49	export SLrep;
50	proc reynolds(poly f){return(SLreynolds(f));};
51	export reynolds;
52	}
53	example
54	{"EXAMPLE:"; echo=2;
55	SL(3);
56	Invar::group;
57	groupideal;
58	print(SLrep);
59	}
60
61	////////////////////////////////////////////////////////////////////////////////
62	// prod(<list>) just gives the product of all entries of <list>.
63	////////////////////////////////////////////////////////////////////////////////
64	static proc prod(list #)
65	"USAGE: prod(<list>)
66	RETURN: the product of all entries of <list>"
67	{ if (size(#)==1) {return(#[1]);};
68	return(#[1]*prod(#[2..size(#)]));
69	};
70
71	////////////////////////////////////////////////////////////////////////////////
72	// monsum(n,d) is the sum of all monomials in x(1),x(2),...,x(n) of degree d
73	////////////////////////////////////////////////////////////////////////////////
74	static proc monsum(int n,int d)
75	"USAGE: monsum(n,d)
76	RETURNS: the sum of all monomials in x(1),x(2),...,x(n) of degree d"
77	{ if (n==1) {return(x(1)^d);};
78	poly output=0;
79	for(int i=0;i<=d;i=i+1)
80	{output=output+x(n)^i*monsum(n-1,d-i);};
81	return(output);
82	};
83
84	////////////////////////////////////////////////////////////////////////////////
85	// sympower(m,d) computes the d-th symmetric power of a representation m
86	////////////////////////////////////////////////////////////////////////////////
87	proc sympower(matrix m,int d)
88	"USAGE: sympower(<matrix>,<int>)
89	RETURNS: If m is a matrix with coefficients in the ring 'group', representing
90	the action on some vector space V, then sympower(m,n) gives the
91	matrix of the representation of the group on the n-th symmetric
92	power of V.
93	EXAMPLE: example sympower; shows an example"
94	{ int n=nrows(m);
95	int l=nvars(group);
96	ring r=0,(z,x(1..n),g(1..l)),(dp(1),dp(n),dp(l));
97	matrix mm=imap(group,m);
98	ideal gideal=std(imap(group,groupideal));
99	matrix vx[n][1]=matrix([x(1..n)]); // vector (x(1),...,x(n))^T
100	poly prodx=prod(x(1..n)); // prodx=x(1)x(2)...*x(n)
101	matrix w[n][1]=mm*vx;
102	map act=r,z,ideal(w);
103	poly ms=monsum(n,d);
104	matrix monlist=coef(ms,prodx); // list of all monomials
105	int j;
106	poly f;
107	matrix q;
108	int s=ncols(monlist); // number of monomials
109	matrix sp[s][s]; // sp : matrix of sym. power
110	for(int i=1;i<=s;i=i+1)
111	{f=monlist[1,i];
112	f=reduce(act(f),gideal)+zms; // zms is added because then
113	q=coef(f,prodx); // all monomials must appear and
114	for(j=1;j<=s;j=j+1) // coef(f,prodx) has right size.
115	{sp[i,j]=q[2,j]-z;}}; // Now subtract the z.
116	setring group;
117	return(imap(r,sp));
118	}
119	example
120	{ "EXAMPLE:"; echo=2;
121	SL(2);
122	print(SLrep);
123	print(sympower(SLrep,3));
124	}
125
126	////////////////////////////////////////////////////////////////////////////////
127	// invar(m) computes the invariant ring of the representation m.
128	////////////////////////////////////////////////////////////////////////////////
129	proc invar(matrix m)
130	"USAGE: invar(<matrix>)
131	RETURNS: If m is a n x n matrix with coefficients in the ring 'group',
132	representing the action on some vector space V, then invar(m);
133	gives polynomials in x(1),x(2),...,x(n) who generate the invariant
134	ring. The following global variables will be set:
135	polyring of type <ring> polynomial ring in x(1),...,x(n)
136	invring of type <ideal> entries generate the inv. ring
137	representation of type <matrix>
138	The base ring will be set to 'polyring' which is a global
139	variable representing the polynomial ring on which the group acts.
140	The variable 'representation' will be set to the input m.
141	EXAMPLE: example invar; shows an example
142	"
143	{
144	int n=nrows(m);
145	int l=nvars(group);
146	matrix representation=m;
147	export representation;
148	ring r=0,(g(1..l),x(1..n),y(1..n)),(dp(l),dp(2*n));
149	matrix mm=imap(group,m);
150	ideal gideal=imap(group,groupideal);
151	matrix vx[n][1]=matrix([x(1..n)]);
152	matrix vy[n][1]=matrix([y(1..n)]);
153	matrix w[1][n]=transpose(vx)*mm-transpose(vy);
154	ideal Gamma=ideal(w),gideal;
155	poly prodg=prod(g(1..l));
156	ideal B=eliminate(Gamma,prodg);
157	print("");
158	print("Ideal B:");
159	print(B);
160	ring polyring=0,(x(1..n)),dp;
161	export polyring;
162	ideal ZFideal=imap(r,B);
163	ZFideal=minbase(ZFideal);
164	print("");
165	print("Zero Fiber Ideal:");
166	print(ZFideal);
167	ideal invring;
168	poly invariant;
169	for(int i=1;i<=size(ZFideal);i=i+1)
170	{invariant=reynolds(ZFideal[i]);
171	if (invariant!=0)
172	{if (invring==0) {invring=invariant;}
173	else {invring=invring,invariant;}}}
174	export invring;
175	print("");
176	print("Generating Invariants:");
177	print(invring);
178	}
179	example
180	{ "EXAMPLE:"; echo=2;
181	SL(2); // Take the group SL_2
182	matrix m=dsum(SLrep,SLrep,SLrep,SLrep);
183	// 4 copies of the standard representation
184	invar(m); // empirical evidence for FFT
185	reynolds(x(1)*x(4)); // The reynolds operator is computed using
186	// the Omega process.
187	}
188
189
190	////////////////////////////////////////////////////////////////////////////////
191	// omega(f,n,i_1,i_2,...,i_t) does the following:
192	// Let M be the matrix of partial derivatives:
193	//
194	// d/d g(1) d/d g(2) ... d/d g(n)
195	// d/d g(n+1) d/d g(n+2) ... d/d g(2n)
196	// . . . .
197	// . . . .
198	// d/d g(n^2-n-1) d/d g(n^2-n+2)... d/d g(n^2)
199	//
200	// Take the submatrix with rows i_1,i_2,...,i_t and columns
201	// 1,2,...,t and apply its determinant (a differential operator) to f.
202	// i_1,i_2,...,i_t is assumed to be descending.
203	//
204	// Cramer's rule is used, which gives us a recursion.
205	////////////////////////////////////////////////////////////////////////////////
206	proc omega(poly f,int n,intvec #)
207	"USAGE: omega(f,n,i_1,i_2,...,i_t)
208	RETURNS: poly (determinant)
209	NOTE:
210	Let M be the matrix of partial derivatives:
211
212	d/d g(1) d/d g(2) ... d/d g(n)
213	d/d g(n+1) d/d g(n+2) ... d/d g(2n)
214	. . . .
215	. . . .
216	d/d g(n^2-n-1) d/d g(n^2-n+2)... d/d g(n^2)
217
218	Take the submatrix with rows i_1,i_2,...,i_t and columns
219	1,2,...,t and apply its determinant (a differential operator) to f.
220	i_1,i_2,...,i_t is assumed to be descending.
221
222	Cramer's rule is used, which gives us a recursion.
223	"
224	{ if (#==0) {return(f);};
225	int m=size(#);
226	if (f==0) {return(0);};
227	if (m==0) {return(f);};
228	poly output=0;
229	intvec a;
230	for(int i=1;i<=m;i=i+1)
231	{if (i==1) {if (m>1) {a=#[i+1..m];};}
232	else {if (m>i) {a=#[1..i-1],#[i+1..m];}
233	else {a=#[1..i-1];};};
234	output=output+(-1)*(i-1)omega(diff(f,g((#[i]-1)*n+m)),n,a);};
235	return(output);
236	};
237
238	////////////////////////////////////////////////////////////////////////////////
239	// SLreynolds(f) applies the reynolds operator to f, if the group is SL_n,
240	// using the omega-process
241	////////////////////////////////////////////////////////////////////////////////
242	proc SLreynolds(poly f)
243	"USAGE: SLreynolds(f) poly f
244	RETURNS: the reynolds operator applied to f
245	NOTE: if the group is SL_n the omega-process is used"
246	{ int nsq=nvars(group);
247	for(int q=1;(q+1)^2<=nsq;q=q+1) {};
248	setring group;
249	int n=nrows(representation);
250	int l=nvars(group);
251	ring r=0,(x(1..n),g(1..l)),(dp(n),dp(l));
252	matrix m=imap(group,representation);
253	ideal gideal=std(imap(group,groupideal));
254	matrix vx[n][1]=matrix([x(1..n)]); // vector (x(1),...,x(n))^T
255	poly prodx=prod(x(1..n)); // prodx=x(1)x(2)...*x(n)
256	matrix w[1][n]=transpose(vx)*m;
257	map act=polyring,ideal(w);
258	poly h=reduce(act(f),gideal);
259	map gtozero=r,x(1..n);
260	poly output=0;number c=1;
261	int i=1;
262	int j;
263	while (h!=0)
264	{output=output+gtozero(h)/c;
265	h=omega(h,q,q..1);
266	for(j=i;j<=i+q-1;j=j+1)
267	{c=c*j;};
268	i=i+1;};
269	setring polyring;
270	return(imap(r,output));
271	};
272
273	////////////////////////////////////////////////////////////////////////////////
274	// torus(n) sets the current group to an n-dimensional torus
275	////////////////////////////////////////////////////////////////////////////////
276	proc torus(int n)
277	"USAGE: torus(<int>)
278	RETURNS: torus(n) sets the current group to an n-dimensional torus. The
279	following global variables will be changed:
280	group of type <ring>
281	groupideal of type <ideal>
282	reynolds of type <proc>
283	The quotient of of `group` and `groupideal` is the coordinate
284	ring of an n-dimensional torus. The basering will be set to
285	'group'.
286	EXAMPLE: example torus; shows an example"
287	{ ring group=0,(g(1..n+1)),dp;
288	export group;
289	ideal groupideal=prod(g(1..n+1))-1;
290	export groupideal;
291	proc reynolds(poly f){return(torusreynolds(f));};
292	export reynolds;
293	}
294	example
295	{"EXAMPLE:"; echo=2;
296	torus(3);
297	Invar::group;
298	groupideal;
299	}
300
301	////////////////////////////////////////////////////////////////////////////////
302	// If m is a list of integer vectors, then torusrep(list m) computes
303	// a matrix with entries in the ring 'group' which represents the
304	// representation of a torus given by the weights m[1],m[2],...
305	////////////////////////////////////////////////////////////////////////////////
306	proc torusrep(list m)
307	"USAGE: torusrep(<list>), <list> must be a list of integer vectors of length
308	n, where n is the dimension of the current torusgroup.
309	RETURNS: torusrep(m) gives a matrix with entries in 'group'. This matrix
310	represents the action of the torus with weights
311	m[1],m[2],...,m[size(m)]
312	EXAMPLE: example torusreynolds; shows an example
313	"
314	{ int r=size(m[1]);
315	int n=size(m);
316	matrix mm[n][n];
317	int min;
318	poly f;
319	int j;
320	for(int i=1;i<=n;i=i+1)
321	{min=0;
322	for(j=1;j<=r;j=j+1)
323	{if (m[i][j]<min) {min=m[i][j];};};
324	f=g(r+1)**(-min);
325	for(j=1;j<=r;j=j+1)
326	{f=fg(j)*(m[i][j]-min);};
327	mm[i,i]=f;};
328	return(mm);
329	}
330	example
331	{ "EXAMPLE:"; echo=2;
332	torus(1); // Take the 1-dimensional torus, the multiplicative group.
333	list weights=-2,-3,7; // 3-dimensial action with weights -2,-3,7
334	matrix m=torusrep(weights);// compute the matrix of the representation
335	invar(m); // compute the invariant ring
336	}
337
338	////////////////////////////////////////////////////////////////////////////////
339	// torusreynolds(f) applies the Reynolds operator to f if the group is
340	// a torus or a finite group.
341	////////////////////////////////////////////////////////////////////////////////
342	proc torusreynolds(poly f)
343	"USAGE: torusreynolds(f) - poly f
344	"RETURNS: the Reynolds operator applied to f if the group is
345	a torus or a finite group.
346	"
347	{ setring group;
348	int n=nrows(representation);
349	int l=nvars(group);
350	ring r=0,(x(1..n),g(1..l)),(dp(n),dp(l));
351	ideal gideal=std(imap(group,groupideal));
352	matrix m=imap(group,representation);
353	matrix vx[n][1]=matrix([x(1..n)]); // vector (x(1),...,x(n))^T
354	matrix w[1][n]=transpose(vx)*m;
355	map act=polyring,ideal(w);
356	poly h=reduce(act(f),gideal);
357	setring polyring;
358	return(imap(r,h));
359	}
360
361	////////////////////////////////////////////////////////////////////////////////
362	// cyclotomic(n,z) gives the n-th cyclotomic polynomial in z, where z is
363	// a variable.
364	////////////////////////////////////////////////////////////////////////////////
365	proc cyclotomic(n,poly z)
366	"USAGE: cyclotomic(n,z) int n, variable z
367	RETURNS: the n-th cyclotomic polynomial in z"
368	{ poly h=z**n-1;
369	ideal id=h;
370	id=std(id);
371	poly f=z-1;
372	int i=prime(n-1);
373	while(i<>2)
374	{if (n%i==0) {f=f(z*(n div i)-1);};
375	i=prime(i-1);};
376	if (n%2==0) {f=f(z*(n div 2)-1);};
377	h=h/gcd(f,h);
378	return(h/leadcoef(h));
379	}
380
381	////////////////////////////////////////////////////////////////////////////////
382	// finite(n) sets the current group to a finite group of order n.
383	////////////////////////////////////////////////////////////////////////////////
384	proc finite(int n)
385	"USAGE: finite(<int>)
386	RETURNS: finite(n) sets the current group to a finite group of order n.
387	The following global variables will be set:
388	group of type <ring>
389	groupideal of type <ideal>
390	reynolds of type <proc>
391	The basering will be set to 'group'.
392	EXAMPLE: example finite; shows an example
393	"
394	{ ring group=0,(g(1),g(2)),dp(2);
395	export group;
396	ideal groupideal=g(2)^n-1,cyclotomic(n,g(1));
397	groupideal=std(groupideal);
398	export groupideal;
399	proc reynolds(poly f){return(torusreynolds(f));};
400	export reynolds; // the procedure torusreynolds
401	// does exactely the right thing
402	// no need to implement
403	// 'finitereynolds'
404	}
405	example
406	{"EXAMPLE:"; echo=2;
407	finite(6); // The symmetric group S_3
408	matrix id=unitmat(3); // identity matrix
409	matrix m3[3][3]=0,1,0,0,0,1,1,0,0; // corresponds with (1 2 3)
410	matrix m2[3][3]=0,1,0,1,0,0,0,0,1; // corresponds with (1 2)
411	list a=id,m3,m3m3,m2,m2m3,m2m3m3; // all elements of S_3
412	matrix rep=finiterep(a); // compute matrix of standard repr.
413	invar(rep); // compute the invariant ring
414	}
415
416
417	////////////////////////////////////////////////////////////////////////////////
418	// If m is a list of matrices, then finiterep(m) gives a matrix with
419	// coefficients in the ring 'group' which represents the action of the
420	// finite group where the elements of the finite group act as
421	// m[1],m[2],...m[size(m)].
422	////////////////////////////////////////////////////////////////////////////////
423	proc finiterep(list m)
424	"USAGE: finiterep(<list>), <list> must be a list of matrices
425	RETURNS: finiterep(m) gives a matrix with coefficients in the ring 'group'
426	which represents the action of the finite group where the elements
427	of the finite group act as m[1],m[2],...m[size(m)].
428	EXAMPLE: example finiterep; shows an example
429	"
430	{ int l=size(m);
431	int n=nrows(m[1]);
432	int i;
433	int j;
434	poly h;
435	matrix finiterep[n][n];
436	for(i=0;i<l;i=i+1)
437	{h=0;
438	for(j=0;j<l;j=j+1)
439	{h=h+g(2)^jg(1)^((ij)%l)/l;};
440	finiterep=finiterep+h*m[i+1];}
441	return(finiterep);
442	}
443	example
444	{"EXAMPLE:"; echo=2;
445	finite(6); // The symmetric group S_3
446	matrix id=unitmat(3); // identity matrix
447	matrix m3[3][3]=0,1,0,0,0,1,1,0,0; // corresponds with (1 2 3)
448	matrix m2[3][3]=0,1,0,1,0,0,0,0,1; // corresponds with (1 2)
449	list a=id,m3,m3m3,m2,m2m3,m2m3m3; // all elements of S_3
450	matrix rep=finiterep(a); // compute matrix of standard repr.
451	invar(rep); // compute the invariant ring
452	}
453

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

Download in other formats:

Original Format