Context Navigation

invar.lib @ 0929198

spielwiese

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

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

Download in other formats:

Original Format