Context Navigation

invar.lib @ f5a8a56

fieker-DuValspielwiese

Last change on this file since f5a8a56 was f5a8a56, checked in by Hans Schoenemann <hannes@…>, 4 years ago
fix: invar.lib: example
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	setring Invar::polyring;
187	Invar::reynolds(x(1)*x(4)); // The reynolds operator is computed using
188	// the Omega process.
189	}
190
191
192	////////////////////////////////////////////////////////////////////////////////
193	// omega(f,n,i_1,i_2,...,i_t) does the following:
194	// Let M be the matrix of partial derivatives:
195	//
196	// d/d g(1) d/d g(2) ... d/d g(n)
197	// d/d g(n+1) d/d g(n+2) ... d/d g(2n)
198	// . . . .
199	// . . . .
200	// d/d g(n^2-n-1) d/d g(n^2-n+2)... d/d g(n^2)
201	//
202	// Take the submatrix with rows i_1,i_2,...,i_t and columns
203	// 1,2,...,t and apply its determinant (a differential operator) to f.
204	// i_1,i_2,...,i_t is assumed to be descending.
205	//
206	// Cramer's rule is used, which gives us a recursion.
207	////////////////////////////////////////////////////////////////////////////////
208	proc omega(poly f,int n,intvec #)
209	"USAGE: omega(f,n,i_1,i_2,...,i_t)
210	RETURNS: poly (determinant)
211	NOTE:
212	Let M be the matrix of partial derivatives:
213
214	d/d g(1) d/d g(2) ... d/d g(n)
215	d/d g(n+1) d/d g(n+2) ... d/d g(2n)
216	. . . .
217	. . . .
218	d/d g(n^2-n-1) d/d g(n^2-n+2)... d/d g(n^2)
219
220	Take the submatrix with rows i_1,i_2,...,i_t and columns
221	1,2,...,t and apply its determinant (a differential operator) to f.
222	i_1,i_2,...,i_t is assumed to be descending.
223
224	Cramer's rule is used, which gives us a recursion.
225	"
226	{ if (#==0) {return(f);};
227	int m=size(#);
228	if (f==0) {return(0);};
229	if (m==0) {return(f);};
230	poly output=0;
231	intvec a;
232	for(int i=1;i<=m;i=i+1)
233	{if (i==1) {if (m>1) {a=#[i+1..m];};}
234	else {if (m>i) {a=#[1..i-1],#[i+1..m];}
235	else {a=#[1..i-1];};};
236	output=output+(-1)*(i-1)omega(diff(f,g((#[i]-1)*n+m)),n,a);};
237	return(output);
238	};
239
240	////////////////////////////////////////////////////////////////////////////////
241	// SLreynolds(f) applies the reynolds operator to f, if the group is SL_n,
242	// using the omega-process
243	////////////////////////////////////////////////////////////////////////////////
244	proc SLreynolds(poly f)
245	"USAGE: SLreynolds(f) poly f
246	RETURNS: the reynolds operator applied to f
247	NOTE: if the group is SL_n the omega-process is used"
248	{ int nsq=nvars(group);
249	for(int q=1;(q+1)^2<=nsq;q=q+1) {};
250	setring group;
251	int n=nrows(representation);
252	int l=nvars(group);
253	ring r=0,(x(1..n),g(1..l)),(dp(n),dp(l));
254	matrix m=imap(group,representation);
255	ideal gideal=std(imap(group,groupideal));
256	matrix vx[n][1]=matrix([x(1..n)]); // vector (x(1),...,x(n))^T
257	poly prodx=prod(x(1..n)); // prodx=x(1)x(2)...*x(n)
258	matrix w[1][n]=transpose(vx)*m;
259	map act=polyring,ideal(w);
260	poly h=reduce(act(f),gideal);
261	map gtozero=r,x(1..n);
262	poly output=0;number c=1;
263	int i=1;
264	int j;
265	while (h!=0)
266	{output=output+gtozero(h)/c;
267	h=omega(h,q,q..1);
268	for(j=i;j<=i+q-1;j=j+1)
269	{c=c*j;};
270	i=i+1;};
271	setring polyring;
272	return(imap(r,output));
273	};
274
275	////////////////////////////////////////////////////////////////////////////////
276	// torus(n) sets the current group to an n-dimensional torus
277	////////////////////////////////////////////////////////////////////////////////
278	proc torus(int n)
279	"USAGE: torus(<int>)
280	RETURNS: torus(n) sets the current group to an n-dimensional torus. The
281	following global variables will be changed:
282	group of type <ring>
283	groupideal of type <ideal>
284	reynolds of type <proc>
285	The quotient of of `group` and `groupideal` is the coordinate
286	ring of an n-dimensional torus. The basering will be set to
287	'group'.
288	EXAMPLE: example torus; shows an example"
289	{ ring group=0,(g(1..n+1)),dp;
290	export group;
291	ideal groupideal=prod(g(1..n+1))-1;
292	export groupideal;
293	proc reynolds(poly f){return(torusreynolds(f));};
294	export reynolds;
295	}
296	example
297	{"EXAMPLE:"; echo=2;
298	torus(3);
299	Invar::group;
300	groupideal;
301	}
302
303	////////////////////////////////////////////////////////////////////////////////
304	// If m is a list of integer vectors, then torusrep(list m) computes
305	// a matrix with entries in the ring 'group' which represents the
306	// representation of a torus given by the weights m[1],m[2],...
307	////////////////////////////////////////////////////////////////////////////////
308	proc torusrep(list m)
309	"USAGE: torusrep(<list>), <list> must be a list of integer vectors of length
310	n, where n is the dimension of the current torusgroup.
311	RETURNS: torusrep(m) gives a matrix with entries in 'group'. This matrix
312	represents the action of the torus with weights
313	m[1],m[2],...,m[size(m)]
314	EXAMPLE: example torusreynolds; shows an example
315	"
316	{ int r=size(m[1]);
317	int n=size(m);
318	matrix mm[n][n];
319	int min;
320	poly f;
321	int j;
322	for(int i=1;i<=n;i=i+1)
323	{min=0;
324	for(j=1;j<=r;j=j+1)
325	{if (m[i][j]<min) {min=m[i][j];};};
326	f=g(r+1)**(-min);
327	for(j=1;j<=r;j=j+1)
328	{f=fg(j)*(m[i][j]-min);};
329	mm[i,i]=f;};
330	return(mm);
331	}
332	example
333	{ "EXAMPLE:"; echo=2;
334	torus(1); // Take the 1-dimensional torus, the multiplicative group.
335	list weights=-2,-3,7; // 3-dimensial action with weights -2,-3,7
336	matrix m=torusrep(weights);// compute the matrix of the representation
337	invar(m); // compute the invariant ring
338	}
339
340	////////////////////////////////////////////////////////////////////////////////
341	// torusreynolds(f) applies the Reynolds operator to f if the group is
342	// a torus or a finite group.
343	////////////////////////////////////////////////////////////////////////////////
344	proc torusreynolds(poly f)
345	"USAGE: torusreynolds(f) - poly f
346	"RETURNS: the Reynolds operator applied to f if the group is
347	a torus or a finite group.
348	"
349	{ setring group;
350	int n=nrows(representation);
351	int l=nvars(group);
352	ring r=0,(x(1..n),g(1..l)),(dp(n),dp(l));
353	ideal gideal=std(imap(group,groupideal));
354	matrix m=imap(group,representation);
355	matrix vx[n][1]=matrix([x(1..n)]); // vector (x(1),...,x(n))^T
356	matrix w[1][n]=transpose(vx)*m;
357	map act=polyring,ideal(w);
358	poly h=reduce(act(f),gideal);
359	setring polyring;
360	return(imap(r,h));
361	}
362
363	////////////////////////////////////////////////////////////////////////////////
364	// cyclotomic(n,z) gives the n-th cyclotomic polynomial in z, where z is
365	// a variable.
366	////////////////////////////////////////////////////////////////////////////////
367	proc cyclotomic(n,poly z)
368	"USAGE: cyclotomic(n,z) int n, variable z
369	RETURNS: the n-th cyclotomic polynomial in z"
370	{ poly h=z**n-1;
371	ideal id=h;
372	id=std(id);
373	poly f=z-1;
374	int i=prime(n-1);
375	while(i<>2)
376	{if (n%i==0) {f=f(z*(n div i)-1);};
377	i=prime(i-1);};
378	if (n%2==0) {f=f(z*(n div 2)-1);};
379	h=h/gcd(f,h);
380	return(h/leadcoef(h));
381	}
382
383	////////////////////////////////////////////////////////////////////////////////
384	// finite(n) sets the current group to a finite group of order n.
385	////////////////////////////////////////////////////////////////////////////////
386	proc finite(int n)
387	"USAGE: finite(<int>)
388	RETURNS: finite(n) sets the current group to a finite group of order n.
389	The following global variables will be set:
390	group of type <ring>
391	groupideal of type <ideal>
392	reynolds of type <proc>
393	The basering will be set to 'group'.
394	EXAMPLE: example finite; shows an example
395	"
396	{ ring group=0,(g(1),g(2)),dp(2);
397	export group;
398	ideal groupideal=g(2)^n-1,cyclotomic(n,g(1));
399	groupideal=std(groupideal);
400	export groupideal;
401	proc reynolds(poly f){return(torusreynolds(f));};
402	export reynolds; // the procedure torusreynolds
403	// does exactely the right thing
404	// no need to implement
405	// 'finitereynolds'
406	}
407	example
408	{"EXAMPLE:"; echo=2;
409	finite(6); // The symmetric group S_3
410	matrix id=unitmat(3); // identity matrix
411	matrix m3[3][3]=0,1,0,0,0,1,1,0,0; // corresponds with (1 2 3)
412	matrix m2[3][3]=0,1,0,1,0,0,0,0,1; // corresponds with (1 2)
413	list a=id,m3,m3m3,m2,m2m3,m2m3m3; // all elements of S_3
414	matrix rep=finiterep(a); // compute matrix of standard repr.
415	invar(rep); // compute the invariant ring
416	}
417
418
419	////////////////////////////////////////////////////////////////////////////////
420	// If m is a list of matrices, then finiterep(m) gives a matrix with
421	// coefficients in the ring 'group' which represents the action of the
422	// finite group where the elements of the finite group act as
423	// m[1],m[2],...m[size(m)].
424	////////////////////////////////////////////////////////////////////////////////
425	proc finiterep(list m)
426	"USAGE: finiterep(<list>), <list> must be a list of matrices
427	RETURNS: finiterep(m) gives a matrix with coefficients in the ring 'group'
428	which represents the action of the finite group where the elements
429	of the finite group act as m[1],m[2],...m[size(m)].
430	EXAMPLE: example finiterep; shows an example
431	"
432	{ int l=size(m);
433	int n=nrows(m[1]);
434	int i;
435	int j;
436	poly h;
437	matrix finiterep[n][n];
438	for(i=0;i<l;i=i+1)
439	{h=0;
440	for(j=0;j<l;j=j+1)
441	{h=h+g(2)^jg(1)^((ij)%l)/l;};
442	finiterep=finiterep+h*m[i+1];}
443	return(finiterep);
444	}
445	example
446	{"EXAMPLE:"; echo=2;
447	finite(6); // The symmetric group S_3
448	matrix id=unitmat(3); // identity matrix
449	matrix m3[3][3]=0,1,0,0,0,1,1,0,0; // corresponds with (1 2 3)
450	matrix m2[3][3]=0,1,0,1,0,0,0,0,1; // corresponds with (1 2)
451	list a=id,m3,m3m3,m2,m2m3,m2m3m3; // all elements of S_3
452	matrix rep=finiterep(a); // compute matrix of standard repr.
453	invar(rep); // compute the invariant ring
454	}
455

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

Download in other formats:

Original Format