1 #include "def.h"
2 #include "macro.h"
3
4 #ifdef DGTRUE /* Darstellungen */
5
6 static INT homtest();
7 static INT tf_idmat();
8 static INT imult();
9 /*
10 static INT nw_nexperm();
11 static INT test_nw_nexperm();
12 static INT nw_nexpart();
13 static INT test_nw_nexpart();
14 */
15
16 /*** maximale Laenge von Permutationen ***/
17 #define TFNMAX 200
18
19 #define FIRST_ST 0
20 #define NEXT_ST 1
21
22
23 /* Zugriff auf das Element (i,j) in der oberen D-Matrix m
24 der Dimension fa; Diagonale ist 1! */
25 #define ODM(m,fa,i,j) m[((long)fa-1L)*(long)(i)-((long)(i)\
26 *((long)(i)+1L))/2L+(j)-1L]
27 /* Zugriff auf das Element (i,j) in der quadratischen Matrix
28 m der Dimension fa */
29 #define MAT(v,fa,i,j) v[(long)fa*(long)(i)+(long)(j)]
30
31
32 /* imult.c */
imult(a,b,d)33 static INT imult(a,b,d) OP a,b,d;
34 /* TF 210989 */
35 /* multipliziert zwei INTEGER matrizen */
36 /* AK 191289 V1.1 */ /* AK 210891 V1.3 */
37 {
38 OP c;
39
40 /* sonderfaelle bei gleichen variablennamen */
41 if (a == d)
42 { c =callocobject(); copy(a,c); imult(c,b,d); freeall(c);
43 return(OK); }
44 if (b == d)
45 { c =callocobject(); copy(b,c); imult(a,c,d); freeall(c);
46 return(OK); }
47
48 /*freigabe des speichers belegt durch d */
49 freeself(d);
50 /*falls beides leere objecte => d auch leer */
51 if (emptyp(a) || emptyp(b)) return(OK);
52
53 if (nullp(a)) return(m_i_i(0L,d));
54 if (nullp(b)) return(m_i_i(0L,d));
55 if (einsp(b)) return(copy(a,d));
56 if (einsp(a)) return(copy(b,d));
57
58 if (S_O_K(a)==MATRIX && S_O_K(b)==MATRIX)
59 return(mult_imatrix_imatrix(a,b,d));
60 else
61 {
62 printobjectkind(a); printobjectkind(b);
63 error("kann ich nicht multiplizieren");
64 return(ERROR);
65 }
66 }
67
68 /**************************************/
69 /*** Spezialfall: zwei INT-Matrizen ***/
70 /**************************************/
mult_imatrix_imatrix(a,b,c)71 INT mult_imatrix_imatrix(a,b,c) OP a,b,c;
72 /* TF 210989 */
73 {
74 #ifdef MATRIXTRUE
75 INT i,j,k,ha,lb,la;
76 OP ll, height,pi,pj;
77 INT zi,aiki,bkji; /* zwischen ergebnis bei matrix-multiplikation */
78
79 if (S_M_LI(a) != S_M_HI(b))
80 {
81 error("matrizen koennen nicht multipliziert werden");
82 return(ERROR);
83 };
84
85 ll=callocobject();
86 height=callocobject();
87 COPY_INTEGER(S_M_H(a),height); ha=S_M_HI(a);
88 COPY_INTEGER(S_M_L(b),ll); lb=S_M_LI(b); la=S_M_LI(a);
89
90 b_lh_m(ll,height,c);
91
92 for (i=0L;i<ha;i++) /* ueber zeilen der linken Matrix */
93 {
94 for (j=0L;j<lb;j++) /* ueber spalten der rechten Matrix */
95 {
96 zi=0L;
97 pi=S_M_IJ(a,i,0L); /* Zeiger auf 1. Objekt in Zeile i */
98 pj=S_M_IJ(b,0L,j); /* Zeiger auf erstes Objekt in Spalte j */
99 for (k=0L; k<la; k++,pi++,pj+=lb)
100 {
101 if ((aiki=S_I_I(pi))==0L) continue;
102 if ((bkji=S_I_I(pj))==0L) continue;
103 zi += aiki*bkji;
104 };
105 m_i_i(zi,S_M_IJ(c,i,j));
106 }
107 }
108 return(OK);
109 #endif
110 }
111
112
test_ndg()113 INT test_ndg()
114 /* zum testen von ndg() und homtest() */
115 /* TF 211189 */
116 /* AK 191289 V1.1 */ /* AK 210891 V1.3 */
117 {
118 INT i,n;
119 OP part = callocobject();
120 OP trans = callocobject();
121 OP nzyk = callocobject();
122 OP transmat = callocobject();
123 OP nzykmat = callocobject();
124 OP ob_n = callocobject();
125
126 scan(PARTITION,part);
127 n=0L;
128 for (i=0L; i<S_PA_LI(part); i++) n+=S_PA_II(part,i);
129 m_i_i(n,ob_n);
130
131 b_ks_p(VECTOR,callocobject(),trans);
132 b_ks_p(VECTOR,callocobject(),nzyk);
133 m_il_v(n,s_p_s(trans));
134 m_il_v(n,s_p_s(nzyk));
135
136 m_i_i(2L,s_p_i(trans,0L));
137 m_i_i(1L,s_p_i(trans,1L));
138 for (i=2L; i < n; i++) m_i_i(i+1L,s_p_i(trans,i));
139 println(trans); ndg(part,trans,transmat); println(transmat);
140
141 for (i=0L; i < n; i++) m_i_i(i+2L,s_p_i(nzyk,i));
142 m_i_i(1L,s_p_i(nzyk,n-1L));
143 println(nzyk); ndg(part,nzyk,nzykmat); println(nzykmat);
144
145 homtest(transmat,nzykmat,ob_n,stdout);
146 freeall(ob_n);freeall(transmat);freeall(nzykmat);freeall(nzyk);
147 freeall(part);freeall(trans);
148 return(OK);
149 }
150
151 /***************************************************************************
152 *
153 * HOMTEST fuer SYMLIB (T. Fuerbringer, August 1989)
154 *
155 * Die Prozedur "homtest" ueberprueft, ob die im Buch von Carmichael
156 * auf Seite 175 (Aufgabe 2) angegebenen Relationen erfuellt sind, d.h.
157 * eine treue Darstellung von Sn erzeugt wird. Es handelt sich insgesamt
158 * um 4+[n/2] Relationen.
159 *
160 * Mit imult linken.
161 *
162 ***************************************************************************/
163
164
165 /***************************************************************************
166 *
167 * Prozedur:
168 * INT tf_idmat(z)
169 *
170 * Beschreibung:
171 * "tf_idmat" testet, ob die Matrix z die Einheitsmatrix ist.
172 *
173 * Rueckgabe:
174 * 1,falls erfuellt oder 0, falls nicht erfuellt.
175 *
176 ***************************************************************************/
tf_idmat(z)177 static INT tf_idmat(z) OP z;
178 /* AK 191289 V1.1 */ /* TF 011289 */ /* AK 210891 V1.3 */
179 {
180 INT i,j,o1,dm;
181
182 dm=S_M_LI(z);
183 for (i=0L; i<dm; ++i)
184 for (j=0L; j<dm; ++j)
185 {
186 o1=S_M_IJI(z,i,j);
187 if (i==j)
188 {
189 if (o1!=1L) return(FALSE);
190 }
191 else
192 {
193 if (o1) return(TRUE);
194 }
195 }
196 return(FALSE);
197 }
198
199 /***************************************************************************
200 *
201 * Prozedur:
202 * INT homtest(trans,nzyk,n,logfp)
203 *
204 * Beschreibung:
205 * Parameter sind
206 * trans : darstellende Matrix einer Transposition
207 * nzyk : darstellende Matrix eines n-Zyklus
208 * n : Sn,
209 * logfp : log-File fuer Zwischenmeldungen (oder NULL)
210 *
211 * Rueckgabe:
212 * Index einer nicht erfuellten Relationen, oder
213 * 0L, falls alle Relationen erfuellt;
214 * Indexreihenfolge:
215 * 1. t^2 = I
216 * 2. s^n = I
217 * 3. (st)^(n-1) = I
218 * 4. (ts^(-1)ts)^3 = I
219 * 3+j. (ts^(-j)ts^j)^2 = I fuer j=2,...,[n/2]
220 *
221 ***************************************************************************/
homtest(trans,nzyk,n,logfp)222 static INT homtest(trans,nzyk,n,logfp) OP trans,nzyk,n; FILE *logfp;
223 /* AK 191289 V1.1 */
224 /* TF 011289 */ /* AK 210891 V1.3 */
225 {
226 OP invzyk = callocobject(),
227 mat = callocobject(),
228 mat1 = callocobject();
229 INT dim,k,i,j,az,ni;
230
231 dim=S_M_LI(trans);
232 ni=S_I_I(n);
233 m_ilih_m(dim,dim,invzyk);
234 m_ilih_m(dim,dim,mat1);
235
236 if (logfp) fprintf(logfp,"%ld Tests:\n",ni/2L+3L);
237 copy(trans,mat);
238 imult(mat,trans,mat1);
239 if (!tf_idmat(mat1))
240 {
241 if (logfp) fprintf(logfp,"Test 1 failed\n");
242 freeall(mat); freeall(mat1); freeall(invzyk);
243 return(1L);
244 }
245 if (logfp)
246 {
247 fprintf(logfp,"Test 1 ok"); fflush(logfp);
248 }
249
250 copy(nzyk,mat); /* (s * t) ** (n-1) = 1 ? */
251 imult(mat,trans,mat);
252 copy(mat,mat1);
253 az=1L;
254 while (2L*az <= (ni-1L))
255 {
256 copy(mat1,invzyk);
257 imult(mat1,invzyk,mat1);
258 az *= 2L;
259 }
260 for (i=az+2L; i<= ni; i++)
261 imult(mat1,mat,mat1);
262 if (!tf_idmat(mat1))
263 {
264 if (logfp) fprintf(logfp,"Test 3 failed\n");
265 freeall(mat); freeall(mat1); freeall(invzyk);
266 return(3L);
267 }
268 if (logfp)
269 {
270 fprintf(logfp,", 3 ok");
271 fflush(logfp);
272 }
273
274 copy(nzyk,mat);
275 az=1L;
276 while (2L*az <= ni-1L)
277 {
278 copy(mat,mat1);
279 imult(mat,mat1,mat);
280 az*=2L;
281 }
282 for (i=az+2L;i<=ni ;++i) /* s**n=1 ? */
283 imult(mat,nzyk,mat);
284 copy(mat,invzyk); /* s**-1 = s**(n-1) */
285 imult(mat,nzyk,mat);
286 if (!tf_idmat(mat))
287 {
288 if (logfp) { fprintf(logfp,"Test 2 failed\n"); fflush(logfp); }
289 freeall(mat); freeall(mat1); freeall(invzyk);
290 return(2L);
291 }
292 if (logfp) { fprintf(logfp,", 2 ok"); fflush(logfp); }
293
294 copy(trans,mat); /* (t*(s**-1)*t*s)**3 = 1 ? */
295 imult(mat,invzyk,mat);
296 imult(mat,trans,mat);
297 imult(mat,nzyk,mat);
298 copy(mat,mat1);
299 imult(mat1,mat,mat1);
300 imult(mat1,mat,mat1);
301 if (!tf_idmat(mat1))
302 {
303 if (logfp) fprintf(logfp,"Test 4 failed \n");
304 freeall(mat); freeall(mat1); freeall(invzyk);
305 return(4L);
306 }
307 if (logfp) { fprintf(logfp,", 4 ok"); fflush(logfp); }
308
309 k=ni/2L; /* ( t * (s**-j) * t * (s**j) ) **2 = 1 fuer j=2,...k ? */
310
311 for (j=2L; j<=k; j++)
312 {
313 imult(mat,nzyk,mat); /* mat ist noch ts^-1 ts */
314 imult(trans,mat,mat);
315 imult(invzyk,mat,mat);
316 imult(trans,mat,mat);
317 copy(mat,mat1);
318 imult(mat1,mat,mat1);
319 if (!tf_idmat(mat1))
320 {
321 if (logfp) fprintf(logfp,"Test %ld failed\n",j+3L);
322 freeall(mat); freeall(mat1); freeall(invzyk);
323 return(j+3L);
324 }
325 if (logfp) { fprintf(logfp,", %ld ok",j+3L); fflush(logfp); }
326 }
327 if (logfp) fprintf(logfp,"\n");
328 freeall(mat); freeall(mat1); freeall(invzyk);
329 return(0L);
330 }
331
332 /*********************************************************************
333 *
334 * Permutationen nach Nijenhuis/Wilf
335 * T. Fuerbringer
336 *
337 *********************************************************************/
338 #ifdef UNDEF
339 /*********************************************************************
340 r=nw_nexperm(m,n,a,sig)
341 INT m : 0 = erste Permutation berechnen
342 1 = folgende Permutation berechnen
343 OP n : Laenge der Permutation
344 OP a : Permutation
345 OP sig: Vorzeichen der berechneten Permutation
346 INT r : 0 = dies war die letzte Permuation
347 1 = noch weitere Permutationen berechenbar
348 *********************************************************************/
nw_nexperm(m,n,a,sig)349 static INT nw_nexperm(m,n,a,sig) INT m; OP n,a,sig;
350 /* AK 191289 V1.1 */
351 /* TF 011289 */ /* AK 210891 V1.3 */
352 {
353 static INT sgn,ni;
354 INT i,j,ai=0,s,d=0;
355
356 if (!m)
357 {
358 ni=S_I_I(n); if (ni<1L) return(0L);
359 b_ks_p(VECTOR,callocobject(),a);
360 m_il_v(ni,s_p_s(a));
361 for (i=0L; i<ni; i++) m_i_i(i+1L,S_P_I(a,i));
362 sgn=1L;
363 if (sig) m_i_i(1L,sig);
364 return((INT)(ni>1L));
365 }
366 if (sig) m_i_i(-sgn,sig);
367 if (sgn<0L)
368 {
369 sgn=1L;
370 s=0L;
371 for (i=1L; i<ni; i++)
372 {
373 ai=S_P_II(a,i);
374 d=0L;
375 for (j=0L; j<i; j++) if (S_P_II(a,j) > ai) d++;
376 s+=d;
377 if (s % 2L && d < i)
378 {
379 for (d=0L; S_P_II(a,d) >= ai; d++);
380 for (j=d+1L; j<i; j++)
381 if (S_P_II(a,j) < ai &&
382 S_P_II(a,j) > S_P_II(a,d)) d=j;
383 break;
384 }
385 else if ((s%2L)==0L && d>0L)
386 {
387 for (d=0L; S_P_II(a,d)<=ai; d++);
388 for (j=d+1L; j<i; j++)
389 if (S_P_II(a,j)>ai &&
390 S_P_II(a,j)<S_P_II(a,d)) d=j;
391 break;
392 }
393 }
394 m_i_i(S_P_II(a,d),S_P_I(a,i));
395 m_i_i(ai,S_P_I(a,d));
396 }
397 else
398 {
399 sgn= -1L;
400 s=S_P_II(a,0L);
401 m_i_i(S_P_II(a,1L),S_P_I(a,0L));
402 m_i_i(s,S_P_I(a,1L));
403 }
404 if (S_P_II(a,ni-1L)!=1L || S_P_II(a,0L)!=(ni%2L+2L)) return(1L);
405 if (ni<=3L) return(0L);
406 j=ni-3L;
407 for (i=0L; i<j; i++) if (S_P_II(a,i+1L)!=S_P_II(a,i)+1L) return(1L);
408 return(0L);
409 }
410
411 #ifdef TEST
412
test_nw_nexperm()413 static INT test_nw_nexperm()
414 /* AK 191289 V1.1 */ /* TF 011289 */ /* AK 210891 V1.3 */
415 {
416 OP a = callocobject();
417 OP n = callocobject();
418 OP s = callocobject();
419 INT e,anz;
420
421 scan(INTEGER,n);
422 e=nw_nexperm(0L,n,a,s);
423 println(a); anz=1L;
424 while (e)
425 {
426 e=nw_nexperm(1L,n,a,s);
427 println(a); anz++;
428 }
429 printf("%ld\n",anz);
430 }
431
432 #endif
433
434
435 /***************************************************************************
436 *
437 * Berechnung der Partitionen zu einem n > 0 nach Nijenhuis/Wilf
438 * T. Fuerbringer
439 *
440 ***************************************************************************/
441
442 /***************************************************************************
443 *
444 * Prozedur:
445 * INT nw_nexpart(mode,n,part)
446 *
447 * Beschreibung:
448 * mode=0:
449 * "nexpart" berechnet erste Partition (n).
450 * Rueckgabe erfolgt in part.
451 *
452 * Statische Variable:
453 * r[0]: Anzahl verschiedener Ziffern in r[1...];
454 * r[i] (i>0): Eine Ziffer der Partition;
455 * m[i] (i>0): Vielfachheit von r[i];
456 * mode<>0:
457 * Die nach Hijenhuis/Wilf auf part folgende Partition wird
458 * berechnet.
459 *
460 * Rueckgabe:
461 * 0: part ist letzte Partition, sonst !=0.
462 *
463 ***************************************************************************/
nw_nexpart(mode,n,part)464 static INT nw_nexpart(mode,n,part) INT mode; OP n,part;
465 /* AK 191289 V1.1 */
466 /* TF 011289 */ /* AK 210891 V1.3 */
467 {
468 INT i,j,d,s,sum,f;
469 static INT ni,r[50],m[50];
470
471 d=r[0];
472 if (mode!=0L)
473 {
474 sum=(r[d]==1L)? m[d--]+1L : 1L;
475 f=r[d]-1L;
476 if (m[d]!=1L)
477 m[d++]--;
478 r[d]=f;
479 m[d]=(sum/f)+1L;
480 s= sum % f;
481 if (s>0L)
482 {
483 r[++d]=s;
484 m[d]=1L;
485 }
486 r[0]=d;
487
488 f=0L;
489 for (i=1L; i<=d; f+=m[i++]);
490
491 m_il_v(f,S_PA_S(part));
492 f=0L;
493 for (i=1L; i<=d; i++)
494 for (j=m[i]; j>0L; j--) m_i_i(r[i],S_PA_I(part,f++));
495 return((INT)(m[d]!=ni));
496 }
497 else
498 {
499 r[0]=m[1]=1L;
500 ni=r[1]=s_i_i(n);
501 b_ks_pa(VECTOR,callocobject(),part);
502 m_il_v(1L,s_pa_s(part));
503 M_I_I(ni,S_PA_I(part,0L));
504 return((INT)(ni!=1L));
505 }
506 }
507
508 #ifdef TEST
509
test_nw_nexpart()510 static INT test_nw_nexpart()
511 /* AK 191289 V1.1 */
512 /* TF 011289 */ /* AK 210891 V1.3 */
513 {
514 OP part = callocobject();
515 OP n = callocobject();
516 INT e;
517
518 scan(INTEGER,n);
519 e=nw_nexpart(0L,n,part);
520 println(part);
521 while(e!=0L)
522 {
523 e=nw_nexpart(1L,n,part);
524 println(part);
525 }
526 }
527
528 #endif
529
530 #endif
531 /************************ 1.TEIL *************************************
532 * bestehend aus PERMFIND.C, NEXTST.C, NATDG.C
533 *********************************************************************/
534
535
536
537 /*** Zeiger auf Bereich fuer Standardtableaux ***/
538 static int *stptr = (int *)0L;
539
540 /*** momentan gueltige Dimension ***/
541 static int stdim = 0;
542
543 /*** Zugriff auf x-tes Tableau ***/
544 #define STAB(x) (stptr+x*(n+1))
545
546 /********************************************************************
547 *
548 * PERMFIND.C
549 *
550 * Prozeduren zur Berechnung von Vertikal- und Horizontalpermutationen
551 * und deren Vorzeichen.
552 *
553 ********************************************************************/
554
555 /********************************************************************
556 *
557 * Prozedur:
558 * void ltmult(pi,t1,t2)
559 *
560 * Beschreibung:
561 * ltmult multipliziert die Permutation pi von links mit
562 * der Transposition (t1 t2).
563 *
564 ********************************************************************/
ltmult(pi,t1,t2)565 static void ltmult( pi , t1 , t2 ) int *pi; int t1; int t2;
566 /* AK 191289 V1.1 */
567 /* TF 011289 */ /* AK 210891 V1.3 */
568 {
569 register int flag = 0; /* flag protokolliert die Anzahl
570 getauschter Elemente */
571
572 while (flag<2)
573 {
574 if (*++pi == t1)
575 { *pi=t2; flag++; }
576 else if (*pi == t2)
577 { *pi=t1; flag++; }
578 }
579 }
580
581 /********************************************************************
582 *
583 * Prozedur:
584 * int sign(pi)
585 *
586 * Beschreibung:
587 * sign berechnet das Vorzeichen der Permutation pi.
588 *
589 * Rueckgabe:
590 * Vorzeichen
591 *
592 ********************************************************************/
sign(pi)593 static int sign(pi ) int *pi;
594 /* AK 191289 V1.1 */
595 /* TF 011289 */ /* AK 210891 V1.3 */
596 {
597 int k,j,i,s,anz;
598 int d[TFNMAX];
599
600 /*** d protokolliert, auf welche Elemente in pi schon
601 zugegriffen wurde ***/
602
603 for (i=1; i<=pi[0]; d[i++]=0);
604
605 /*** Vorzeichen berechnen ***/
606
607 s=1; /* sign := 1 */
608 k=1; /* Elementindex, ab dem gesucht wird */
609 anz=0; /* Anzahl erledigter Elemente */
610 while (anz < pi[0])
611 {
612 i=k++;
613 while (d[i]) i++; /* suche unberuehrtes Element */
614 j=i;
615 d[i]=1;
616 anz++;
617 while (pi[i]!=j) /* durchlaufe zu i gehoerenden Zykel */
618 {
619 s= -s; /* Vorzeichen dreht sich staendig um */
620 i=pi[i];
621 d[i]=1; /* alle Zykelelemente als erledigt kennzeichnen */
622 anz++;
623 if (k==i) k++; /* vor i muss nicht mehr gesucht werden */
624 }
625 }
626 return(s);
627 }
628
629 /********************************************************************
630 *
631 * Prozedur:
632 * int find_pq(alpha,t1,t2,q,ps)
633 *
634 * Beschreibung:
635 * Seien t1 und t2 Tableaux zum Rahmen [alpha].
636 * find_pq sucht q aus V(t1), ps aus H(q*t1) mit t2= ps*q t1.
637 *
638 * Rueckgabe:
639 * 0 gefunden, -1 sonst.
640 *
641 ********************************************************************/
find_pq(alpha,t1,t2,q,ps)642 static int find_pq(alpha,t1,t2,q,ps) int *alpha,*t1,*t2; int *q; int *ps;
643 /* AK 191289 V1.1 */
644 /* TF 011289 */ /* AK 210891 V1.3 */
645 {
646 register int ss,zs,ks,n,i;
647 int tloc[TFNMAX];
648 int d[TFNMAX],z,k,ks_pos,k_pos;
649 int zpos[TFNMAX],spos[TFNMAX];
650
651 n=tloc[0]=ps[0]=q[0]=t1[0];
652
653 /*** erstmal vorbesetzen ***/
654
655 for (i=1; i<=n; i++)
656 {
657 d[i]=0; /* Zugriffsprotokoll auf Elemente in t2 */
658 q[i]=i; /* q := id */
659 tloc[i]=t1[i]; /* tloc:=t1, da in t1 selbst keine Elemente
660 getauscht werden duerfen */
661 }
662
663 /*** bestimme Zeilen- und Spaltenpositionen
664 der Elemente im Tableau tloc ***/
665
666 i=1; /* Laufindex durchs Tableau tloc */
667 for (zs=1; zs<=alpha[0]; zs++) /* Zeilen */
668 for (ss=1; ss<=alpha[zs]; ss++) /* Spalten */
669 {
670 ks=tloc[i++];
671 zpos[ks]=zs;
672 spos[ks]=ss;
673 }
674
675 /*** durchlaufe nun t2 v.l.n.r. und v.o.n.u.: Position (zs,ss) ***/
676
677 ks_pos=1; /* Position des Laufelements im Vektor t2 */
678 k_pos=0; /* Position von Element k in tloc */
679 for (zs=1; zs<=alpha[0]; zs++)
680 {
681 for (ss=1; ss<=alpha[zs]; ss++)
682 {
683 ks=t2[ks_pos++]; /* ks ist das Element in Zeile zs, Spalte ss */
684 z=zpos[ks]; /* z ist Zeilenposition von ks in tloc */
685 if (z != zs) /* z!=zs, dann muss man in tloc vertauschen */
686 {
687 k=tloc[k_pos+spos[ks]]; /* Vertauschungselement in tloc an
688 Position (zs,s), wobei s die Spalten-
689 position von ks in tloc ist */
690 if (d[k]) return(-1); /* Fehler: Element wurde schon mal
691 getauscht: kein q gefunden */
692 ltmult(q,k,ks); /* sonst: q = (k ks)q */
693 ltmult(tloc,k,ks); /* und tloc ebenso */
694 zpos[ks]=zs; /* auch die Zeilenpositionen sind neu */
695 zpos[k]=z;
696 }
697 d[ks]=1; /* ks erledigt */
698 }
699 k_pos+=alpha[zs]; /* k_pos zeigt jetzt auf 1. Element in der
700 naechsten Zeile von tloc */
701 }
702
703 /*** Nun muss ggf. noch ps bestimmt werden. Das ist aber fuer
704 dieses Programm nicht notwendig. Fuer andere Anwendungen
705 sind die Kommentarklammern zu entfernen ***/
706
707 /* for (i=1; i<=n; i++) ps[tloc[i]]=t2[i]; */
708
709 return(0);
710 }
711
712 /********************************************************************
713 *
714 * NEXTST.C
715 *
716 * Prozeduren zur Berechnung der Standardtableaux
717 * (geordnet nach der letzten Ziffer)
718 *
719 ********************************************************************/
720
721 /********************************************************************
722 *
723 * Prozedur:
724 * void set_LL_min(st,talpha,alpha)
725 *
726 * Beschreibung:
727 * set_LL_min belegt das Tableau zum Teilrahmen talpha von alpha
728 * so vor, dass es im Teilrahmen LL-kleinstes ist. st ist das zu
729 * alpha gehoerende Tableau.
730 *
731 ********************************************************************/
set_LL_min(st,talpha,alpha)732 static void set_LL_min(st,talpha,alpha) int *st; int *talpha; int *alpha;
733 /* AK 191289 V1.1 */
734 /* TF 011289 */ /* AK 210891 V1.3 */
735 {
736 register int k,s,z,*pos;
737
738 /*** k ist die einzusetzende Ziffer, d.h k=1,2,... ***/
739
740 k=1;
741
742 /*** Durchlaufe nun den Teilrahmen spalten- und zeilenweise ***/
743
744 for (s=1; s<=talpha[1]; s++)
745 {
746 pos=st+s; /* "absolute" (Vektor-)Position in st */
747 for (z=1; talpha[z] >= s; z++)
748 {
749 *pos=k++;
750 pos+=alpha[z]; /* eine Zeile weiter */
751 if (z >= talpha[0]) break;
752 }
753 }
754 }
755
756 /********************************************************************
757 *
758 * Prozedur:
759 * int nextst(mode,alpha,st)
760 *
761 * Beschreibung:
762 *
763 * mode=FIRST_ST:
764 * nextst berechnet das LL-kleinste Standardtableau.
765 *
766 * mode=NEXT_ST:
767 * nextst berechnet das auf st folgende Tableau in der LLS.
768 *
769 * Rueckgabe:
770 * 0 zeigt an, dass st schon letztes war, sonst 1.
771 *
772 ********************************************************************/
nextst(mode,alpha,st)773 static int nextst(mode,alpha,st) int mode; int *alpha; int *st;
774 /* AK 191289 V1.1 */
775 /* TF 011289 */ /* AK 210891 V1.3 */
776 {
777 int talpha[TFNMAX];
778 register int n,s,z,zz,pos,i,j,k;
779
780 /*** LL-kleinstes ***/
781
782 if (mode==FIRST_ST)
783 {
784 st[0]=0; /* n bestimmen */
785 for (i=1; i<=alpha[0]; i++) st[0]+=alpha[i];
786 set_LL_min(st,alpha,alpha);
787 return(1);
788 }
789 else
790
791 /*** LL-naechstes ***/
792
793 {
794
795 /*** berechne kleinstes Teiltableau, so dass man die
796 groesste Ziffer darin nach unten schieben kann ***/
797
798 n=st[0];
799 talpha[0]=1;
800 talpha[1]=1; /* Teiltableau talpha mit [1] vorbesetzen */
801 for (i=2; i<=n; i++) /* durchlaufe alle Ziffern 2,...,n */
802 {
803
804 /*** bestimme Zeile z und "absolute" Position k von i in st ***/
805
806 z=1; /* (1,1) ist immer mit 1 besetzt -> uebergehen */
807 s=2;
808 k=1;
809 while (st[++k]!=i)
810 {
811 if (s++>=alpha[z])
812 {
813 s=1;
814 z++;
815 }
816 }
817
818 /*** erweitere Rahmen talpha um Kaestchen fuer i ***/
819
820 if (z > talpha[0])
821 {
822 talpha[0]++;
823 talpha[z]=0;
824 }
825 talpha[z]++;
826
827 /*** pruefe nun in talpha, ob man i nach unten bringen kann ***/
828
829 if (talpha[0] > 1 && talpha[1] > 1) /* talpha gross genug ? */
830 {
831 pos=k+alpha[z]-s;
832 for (zz=z+1; zz <= talpha[0]; zz++)
833 /* die Zeilen unter z testen */
834 {
835 pos+=talpha[zz]; /* pos zeigt auf Endziffer in Zeile zz */
836 if (zz==talpha[0]) /* Ist zz die letzte Zeile ... */
837 j=1;
838 else if (talpha[zz] > talpha[zz+1])
839 /* oder steht nichts darunter, ... */
840 j=1; /* so ist tauschen mglich */
841 else
842 j=0; /* sonst nicht */
843 if (j && st[pos] < i) /* Tauschen erlaubt, falls Endziffer
844 kleiner als Ziffer i ist */
845 {
846 s=st[pos]; /* tausche i mit st[pos] */
847 st[pos]=i;
848 st[k]=s;
849 if (--talpha[zz]==0) talpha[0]--;
850 /* Teiltableau ohne Ziffer i neu sortieren */
851 if (talpha[0] > 1 && talpha[1] > 1)
852 set_LL_min(st,talpha,alpha);
853 return(1); /* fertig */
854 }
855 pos+=alpha[zz]-talpha[zz];
856 /* sonst naechste Zeile testen */
857 }
858 }
859 }
860
861 /*** alle Ziffern getestet und keine Tauschmoeglichkeit gefunden:
862 letztes Tableau erreicht ***/
863
864 return(0);
865 }
866 }
867
868 /********************************************************************
869 *
870 * NATDG.C (ohne getdim)
871 *
872 * Prozeduren zur Berechnung der natuerlichen Darstellung der Sn
873 * Version mit Pufferung aller Standardtableaux.
874 *
875 ********************************************************************/
876
877 /********************************************************************
878 *
879 * Prozedur:
880 * int allst(alpha,n)
881 *
882 * Beschreibung:
883 * allst berechnet alle Standardtableaux, belegt entsprechend
884 * viel Speicher und liefert Zeiger auf den Speicherbereich.
885 * Der Bereich muss mit free(Bereich) freigegeben werden.
886 *
887 * Rueckgabe:
888 * 0 (OK) oder <>0 (FEHLER)
889 ********************************************************************/
allst(alpha,n)890 static int allst(alpha,n) int *alpha; int n;
891 /* AK 191289 V1.1 */
892 /* TF 011289 */ /* AK 210891 V1.3 */
893 {
894 int i,e,*b;
895 int st[TFNMAX];
896
897 if (stptr)
898 {
899 SYM_free(stptr);
900 stptr=(int *)0L;
901 }
902 stptr=b=(int *)SYM_malloc(sizeof(int)*(long)(n+1)*(long)stdim);
903 if (!b) return(-1);
904 e=nextst(FIRST_ST,alpha,st);
905 while (e)
906 {
907 for (i=0; i<=n; *b++=st[i++]);
908 e=nextst(NEXT_ST,alpha,st);
909 }
910 return(0);
911 }
912
913 /********************************************************************
914 *
915 * Prozedur:
916 * int koeff(alpha,pi,stk,stj)
917 *
918 * Beschreibung:
919 * koeff berechnet den Koeffizienten der Permutation pi
920 * im Gruppenalgebra-Element e(kj) := e(k)*pi(kj).
921 * Dabei beziehen sich k und j auf die zugehoerigen LLS-geordneten
922 * Standardtableaux stk und stj, wobei stj = pi(kj)stk.
923 *
924 * Rueckgabe:
925 * Koeffizient -1,0,1
926 *
927 ********************************************************************/
koeff(alpha,pi,stk,stj)928 static int koeff(alpha,pi,stk,stj) int *alpha; int *pi; int *stk; int *stj;
929 /* AK 191289 V1.1 */
930 /* TF 011289 */
931 {
932 int i,ps[TFNMAX],q[TFNMAX];
933 int t2[TFNMAX];
934
935 /*** berechne Tableau t2 = pi * st[j] ***/
936
937 t2[0]=pi[0];
938 for (i=1; i<=pi[0]; i++)
939 t2[i]=pi[stj[i]];
940
941 /*** finde nun Permutationen q,ps ***/
942
943 if (find_pq(alpha,stk,t2,q,ps)<0) return(0);
944
945 /*** der Koeffizient ist dann sgn(q) ***/
946
947 return(sign(q));
948 }
949
950 /********************************************************************
951 *
952 * Prozedur:
953 * int koeffid(alpha,stk,stj)
954 *
955 * Beschreibung:
956 * koeffid ist ein Spezialfall von koeff fuer pi=id.
957 *
958 * Rueckgabe:
959 * Koeffizient -1,0,1
960 *
961 ********************************************************************/
koeffid(alpha,stk,stj)962 static int koeffid(alpha,stk,stj) int *alpha; int *stk; int *stj;
963 /* AK 191289 V1.1 */
964 /* TF 011289 */ /* AK 210891 V1.3 */
965 {
966 int ps[TFNMAX],q[TFNMAX];
967
968 /*** finde nun Permutationen q,ps ***/
969
970 if (find_pq(alpha,stk,stj,q,ps) < 0) return(0);
971
972 /*** der Koeffizient ist dann sgn(q) ***/
973
974 return(sign(q));
975 }
976
977 /********************************************************************
978 *
979 * Prozedur:
980 * int ndg_L_alpha(alpha,la)
981 *
982 * Beschreibung:
983 * ndg_L_alpha berechnet die Matrix L_alpha := M_alpha ^ (-1),
984 * mit M_alpha = (( koeff(id,e(kj) )).
985 * Die Matrix L_alpha wird in la zurueckgegeben.
986 * Es wird nur die rechte obere Haelfte der oberen Dreiecksmatrix
987 * gespeichert. Auch die Diagonale aus Einsen wird nicht gespeichert.
988 * Die Matrix liegt dann in Vektorform vor. Auf die Komponenten
989 * kann mit dem Makro ODM(la,z,s) zugegriffen werden.
990 *
991 * Rueckgabe:
992 * 0, -1 fuer Fehler
993 *
994 ********************************************************************/
ndg_L_alpha(alpha,la)995 static int ndg_L_alpha(alpha,la) int *alpha; int *la;
996 /* AK 191289 V1.1 */
997 /* TF 011289 */ /* AK 210891 V1.3 */
998 {
999 int dim,n,*TF_la,*TF_la1,*TF_la2,o1,o2,i,k;
1000 int id[TFNMAX];
1001 int j;
1002 int stj,stk;
1003
1004 /*** erstmal n aus der Partition berechnen ***/
1005
1006 n=0;
1007 for (i=1; i<=alpha[0]; n+=alpha[i++]);
1008
1009 id[0]=n;
1010 for (i=1; i<=n; i++) id[i]=i;
1011
1012 TF_la=la;
1013
1014 /*** alle Tableaux berechnen, dim muss bekannt sein ***/
1015
1016 if (stptr)
1017 {
1018 SYM_free(stptr);
1019 stptr=(int *)0L;
1020 }
1021 dim=stdim;
1022 if (allst(alpha,n)) return(-1);
1023
1024 /*** M_alpha transponiert berechnen ***/
1025
1026 for (stj=0; stj < dim ; stj++)
1027 for (stk=stj+1; stk < dim; stk++)
1028 *TF_la++=koeffid(alpha,STAB(stk),STAB(stj));
1029
1030 /**** dann L_alpha := M_alpha hoch -1 berechnen ***/
1031
1032 for (j=dim-1; j>=1; j--)
1033 {
1034 TF_la=la+dim*(long)j-((long)j*((long)j+1L))/2L-(long)j-1L;
1035 TF_la1=TF_la2= &ODM(la,dim,j-1,j);
1036 for (i=j-1; i>=0; i--)
1037 {
1038 *TF_la1= -(*TF_la1); /* ODM(la,dim,i,j) */
1039 TF_la1 -= dim-i-1; /* eine Zeile aufwaerts */
1040 }
1041 for (k=j+1; k<dim; k++)
1042 {
1043 if ((o2=TF_la[k])==0) continue; /* ODM(la,dim,j,k)==0 ? */
1044 TF_la1=TF_la2; /* zeigt auf ODM(la,dim,j-1,j) */
1045 for (i=j-1; i >= 0; i--)
1046 {
1047 if ((o1= *TF_la1)!=0) /* ODM(la,dim,i,j) */
1048 ODM(la,dim,i,k) += o1*o2;
1049 TF_la1 -= dim-i-1; /* eine Zeile hoch */
1050 }
1051 }
1052 }
1053 return(0);
1054 }
1055
1056 /********************************************************************
1057 *
1058 * Prozedur:
1059 * int ndg_P_pi(alpha,pi,p)
1060 *
1061 * Beschreibung:
1062 * ndg_P_pi berechnet die Koeffizientenmatrix P_pi zur Permutation
1063 * pi.
1064 * Die Matrix wird wie bei L_alpha als Vektor gespeichert.
1065 * Zugriff erfolgt mit dem Makro MAT(pi,z,s).
1066 *
1067 * Rueckgabe:
1068 * 0
1069 *
1070 ********************************************************************/
ndg_P_pi(alpha,pi,p)1071 static int ndg_P_pi(alpha,pi,p) int *alpha; int *pi; int *p;
1072 /* AK 191289 V1.1 */
1073 /* TF 011289 */ /* AK 210891 V1.3 */
1074 {
1075 int n,k,*TF_p;
1076 int pi1[TFNMAX];
1077 int stj,stk;
1078
1079 TF_p=p;
1080
1081 /*** berechne pi1 = pi hoch -1 ***/
1082
1083 n=pi1[0]=pi[0];
1084 for (k=1; k<=n; k++) pi1[pi[k]]=k;
1085
1086 /*** berechne Koeffizientenmatrix ***/
1087
1088 for (stk=0; stk < stdim; stk++)
1089 for (stj=0; stj < stdim; stj++)
1090 *TF_p++=koeff(alpha,pi1,STAB(stj),STAB(stk));
1091
1092 return(0);
1093 }
1094
1095 /********************************************************************
1096 *
1097 * Prozedur:
1098 * int tfmult(la,p,dim,fp,dig)
1099 *
1100 * Beschreibung:
1101 * tfmult berechnet Produkt p=la*p. Dabei wird die Dreiecksstruktur
1102 * bei la zugrunde gelegt, ebenso die Vektorform der Matrizen.
1103 * Als Ergebnis erhaelt man die Darstellungsmatrix.
1104 *
1105 * Ist fp!=NULL, so werden Werte mit Breite <dig> auf FILE
1106 * fp ausgegeben. Bei dig=0 werden die Werte nur durch
1107 * ein Leerzeichen getrennt.
1108 *
1109 * Rueckgabe:
1110 * Charakterwert fuer P(pi)
1111 *
1112 ********************************************************************/
tfmult(la,p,dim,fp,dig)1113 static int tfmult(la,p,dim,fp,dig) int *la; int *p; int dim; FILE *fp; int dig;
1114 /* AK 191289 V1.1 */
1115 /* TF 011289 */ /* AK 210891 V1.3 */
1116 {
1117 int o1,p2,*TF_la,*TF_p,*TF_pp,i,j;
1118 int h,k,ch;
1119 char format[6];
1120
1121 /*** Format-String fuer die Ausgabe ***/
1122
1123 if (dig) sprintf(format,"%%%dd",dig);
1124 else strcpy(format,"%d ");
1125
1126 /*** multipliziere Matrizen ***/
1127
1128 ch=0; /* Charakter = 0 */
1129 for (j=0; j<dim; j++)
1130 {
1131 TF_la=la;
1132 TF_p=p+(long)j;
1133 for (i=0; i<dim; i++)
1134 {
1135 h= *TF_p; /* MAT(p,dim,i,j) */
1136 TF_pp=TF_p+(long)dim;
1137 for (k=i+1; k<dim; k++,TF_pp+=(long)dim)
1138 {
1139 if ((o1= *TF_la++)==0) continue; /* ODM(i,k) */
1140 if ((p2= *TF_pp) != 0) /* MAT(p,dim,k,j)*/
1141 h = (p2 > 0) ? h+o1 : h-o1 ;
1142 }
1143 *TF_p=h; /* MAT(p,dim,i,j) */
1144 if (i==j) ch+=h;
1145 TF_p+=dim;
1146 }
1147 }
1148
1149 /*** Ausgabe ***/
1150
1151 if (fp!=NULL)
1152 {
1153 TF_la=p;
1154 for (i=0; i<dim; i++)
1155 {
1156 for (j=0; j<dim; j++) fprintf(fp,format,*TF_la++);
1157 fprintf(fp,"\n");
1158 }
1159 }
1160
1161 return(ch);
1162 }
1163
1164 #ifdef UNDEF
1165 /********************************************************************
1166 *
1167 * Prozedur:
1168 * int tfchmult(alpha,la,pi,dim)
1169 *
1170 * Beschreibung:
1171 * tfchmult berechnet Charakter mittels dem Produkt la*p(pi).
1172 * Die Matrix p(pi) wird berechnet, aber nicht gespeichert.
1173 *
1174 * Rueckgabe:
1175 * Charakterwert
1176 *
1177 ********************************************************************/
tfchmult(alpha,la,pi)1178 static int tfchmult(alpha,la,pi) int *la; int *pi; int *alpha;
1179 /* AK 191289 V1.1 */
1180 /* TF 011289 */ /* AK 210891 V1.3 */
1181 {
1182 int h,ch,k,n,o1,p2;
1183 int *TF_la,pi1[TFNMAX],stj,stk;
1184
1185 /*** berechne Charakter ***/
1186
1187 ch=0;
1188 TF_la=la;
1189
1190 /*** berechne pi1 = pi hoch -1 ***/
1191
1192 n=pi1[0]=pi[0];
1193 for (k=1; k<=n; k++) pi1[pi[k]]=k;
1194
1195 /*** durchlaufe alle Zeilen von la ***/
1196
1197 for (stk=0; stk < stdim; stk++)
1198 {
1199
1200 /*** durchlaufe entsprechende Spalte in p(pi) ***/
1201
1202 h=koeff(alpha,pi1,STAB(stk),STAB(stk)); /* Diagonalelement */
1203 for (stj=stk+1; stj < stdim; stj++)
1204 {
1205 if ((o1= *TF_la++)!=0)
1206 {
1207 p2=koeff(alpha,pi1,STAB(stk),STAB(stj));
1208 if (p2) h+=o1*p2;
1209 }
1210 }
1211 ch+=h;
1212 }
1213 return(ch);
1214 }
1215 #endif
1216 /****************************** TEIL 2 *******************************
1217 * hier befindet sich die Bindung zu SYMCHAR
1218 *********************************************************************/
1219
1220 /*********************************************************************
1221 *
1222 * Prozedur:
1223 * ndg(part,perm,dg_mat)
1224 *
1225 * Beschreibung:
1226 * "ndg" berechnet Darstellungs-Matrix dg_mat
1227 * zur Permutation perm und Partition part.
1228 *
1229 ****************************************************************************/
ndg(part,perm,dg_mat)1230 INT ndg(part,perm,dg_mat) OP part,perm,dg_mat;
1231 /* AK 191289 V1.1 */ /* TF 011289 */ /* AK 210891 V1.3 */
1232 {
1233 int alpha[TFNMAX],pi[TFNMAX],*la,*p,*pp;
1234 INT i,j;
1235 long dim;
1236 OP opdim = callocobject();
1237
1238 dimension(part,opdim);
1239 dim=(long)S_I_I(opdim);
1240 stdim=(int)dim;
1241 freeall(opdim);
1242
1243 la=(int *)SYM_calloc((dim*(dim-1L))/2L + 1L,sizeof(int));
1244 p=(int *)SYM_calloc(dim*dim,sizeof(int));
1245
1246 if (!p)
1247 {
1248 error("no memory available for matrices");
1249 return(ERROR);
1250 }
1251 if ((!la) && (dim > 1L))
1252 {
1253 error("no memory available for matrices");
1254 return(ERROR);
1255 }
1256 alpha[0]=(int)S_PA_LI(part);
1257 for (i=0L; i<S_PA_LI(part); i++)
1258 alpha[alpha[0]-i]=(int)S_PA_II(part,i);
1259 pi[0]=(int)S_P_LI(perm);
1260 for (i=0L; i<S_P_LI(perm); i++)
1261 pi[i+1]=(int)S_P_II(perm,i);
1262
1263 if (ndg_L_alpha(alpha,la))
1264 {
1265 error("calling ndg_L_alpha / allst");
1266 return(ERROR);
1267 }
1268 ndg_P_pi(alpha,pi,p);
1269 tfmult(la,p,(int)dim,NULL,0);
1270 SYM_free(la);
1271 if (stptr)
1272 {
1273 SYM_free(stptr);
1274 stptr=(int *)0L;
1275 }
1276 pp=p;
1277 m_ilih_m((INT)dim,(INT)dim,dg_mat);
1278 for (i=0L; i<dim; i++)
1279 for (j=0L; j<dim; j++) m_i_i((INT)(*pp++),S_M_IJ(dg_mat,i,j));
1280 SYM_free(p);
1281 return(OK);
1282 }
1283
1284 #endif /* DGTRUE */
1285
1286