1 /* IK: 150591 */
2 #include "def.h"
3 #include "macro.h"
4 static INT berechne_Dominanz();
5 static INT berechne_sprod_kk();
6 static INT berechne_sumvec();
7 static INT berechne_M_Pk();
8 static INT mult_udrmatrix_matrix();
9 
10 #ifdef MATRIXTRUE
11 #ifdef CHARTRUE
compute_zonal_with_alphabet(part,l,res)12 INT compute_zonal_with_alphabet(part,l,res) OP part,l,res;
13 /* AK 210891 V1.3 */
14 {
15 	OP c,d,e,f,g,h,i;
16 	INT ind,z;
17 	INT erg = OK;
18 
19 	C2R(part,l,"compute_zonal_with_alphabet",res);
20 	CTO(PARTITION,"compute_zonal_with_alphabet",part);
21 	CTO(INTEGER,"compute_zonal_with_alphabet",l);
22 	if (S_PA_LI(part) > S_I_I(l))
23 		{
24 		erg += init(POLYNOM,res);
25 		goto s2r_ende;
26 		}
27 	c = callocobject();
28 	d = callocobject();
29 	e = callocobject();
30 	f = callocobject();
31 	g = callocobject();
32 	h = callocobject();
33 	i = callocobject();
34 
35 	erg += weight(part,c);
36 	erg += makevectorofpart(c,e);
37 
38 	erg += young_tafel(c,h,NULL,NULL);
39 	erg += invers(h,h);
40 	erg += transpose(h,i);
41 
42 	erg += m_ilih_m(S_V_LI(e),S_V_LI(e),f);
43 	erg += berechne_Dominanz(S_V_LI(e),e,f);
44 
45 	erg += m_ilih_m(S_V_LI(e),S_V_LI(e),g);
46 	erg += berechne_sprod_kk(S_V_LI(e),c,e,h,i, f, g);
47 
48 	erg += m_ilih_m(S_V_LI(e),S_V_LI(e),d);
49 	erg += berechne_M_Pk(S_V_LI(e),g,e,f,d);
50 
51 	ind = indexofpart(part);
52 	erg += init(POLYNOM,res);
53 	for (z=0L;z<S_V_LI(e);z++)
54 		{
55 		if (not nullp(S_M_IJ(d,ind,z)))
56 			{
57 			erg += compute_monomial_with_alphabet(S_V_I(e,z),l,f);
58 			erg += mult_apply(S_M_IJ(d,ind,z),f);
59 			erg += add_apply(f,res);
60 			}
61 		}
62 
63 	erg += freeall(c);
64 	erg += freeall(d);
65 	erg += freeall(e);
66 	erg += freeall(f);
67 	erg += freeall(g);
68 	erg += freeall(h);
69 	erg += freeall(i);
70 s2r_ende:
71 	S2R(part,l,"compute_zonal_with_alphabet",res);
72 	ENDR("compute_zonal_with_alphabet");
73 }
74 #endif /* CHARTRUE */
75 #endif /* MATRIXTRUE */
76 
77 
78 
79 /*-----------------------------------------------------------------*/
80 /*                         UNTERPROGRAMME                          */
81 /*-----------------------------------------------------------------*/
82 
83 /*-----------------------------------------------------------------*/
84 #ifdef PARTTRUE
berechne_sumvec(dim,vectorofpart,laenge,sumvec)85 static INT berechne_sumvec(dim,vectorofpart,laenge,sumvec) INT dim;
86 	OP vectorofpart,laenge,sumvec;
87               /* sumvec enthaelt an der Stelle i den zu der Partition
88                  mit dem Index i gehoerenden Teilsummenvektor, also
89                  den Vektor der als i-ten Eintrag die Summe der
90                  ersten i Teile der Partition enthaelt */
91 
92 /* AK 210891 V1.3 */
93 {
94 	INT i,j,k;
95 
96 	for (i=1L;i<dim-1L;i++)
97 	   {
98 	   m_l_v(S_V_I(laenge,i),S_V_I(sumvec,i));
99 	   for (j=0L;j<S_V_II(laenge,i);j++)
100 	      {
101 	      m_i_i(0L,S_V_I(S_V_I(sumvec,i),j));
102 	      for (k=0L;k<=j;k++)
103 		 {
104 	add(S_V_I(S_V_I(sumvec,i),j),S_PA_I(S_V_I(vectorofpart,i),
105 	    S_V_II(laenge,i)-1-k),S_V_I(S_V_I(sumvec,i),j));
106 		 }
107 	      } /* for j */
108 	   } /* for i */
109 	return OK;
110 }                                          /* Ende berechne_sumvec */
111 #endif /* PARTTRUE */
112 /*-----------------------------------------------------------------*/
113 
114 
115 
116 
117 /*-----------------------------------------------------------------*/
118 /*-----------------------------------------------------------------*/
119 /*             siehe Programmbeschreibung unter 2.                 */
120 /*-----------------------------------------------------------------*/
121 #ifdef MATRIXTRUE
122 #ifdef PARTTRUE
berechne_Dominanz(dim,vectorofpart,Dominanz)123 static INT berechne_Dominanz(dim,vectorofpart,Dominanz) INT dim;
124 	OP vectorofpart,Dominanz;
125 /* AK 210891 V1.3 */
126 {
127 	INT i,j,k;
128 	OP laenge,sumvec;
129 
130 	laenge=callocobject();
131 	sumvec=callocobject();
132 
133 	for (i=1L;i<dim-1L;i++)
134 	   {
135 	   m_i_i(1L,S_M_IJ(Dominanz,i,i));
136 	   m_i_i(1L,S_M_IJ(Dominanz,0L,i));
137 	   m_i_i(1L,S_M_IJ(Dominanz,i,dim-1L));
138 	   }
139 	m_i_i(1L,S_M_IJ(Dominanz,0L,0L));
140 	m_i_i(1L,S_M_IJ(Dominanz,0L,dim-1L));
141 	m_i_i(1L,S_M_IJ(Dominanz,dim-1L,dim-1L));
142 	m_il_v(dim,laenge);       /* laenge enthaelt an Stelle i die Laenge
143 				     der zum Index i gehoerenden Partition */
144 
145 	for (i=0L;i<dim;i++)
146 	   M_I_I(S_PA_LI(S_V_I(vectorofpart,i)),S_V_I(laenge,i));
147 	m_il_v(dim,sumvec);
148 	berechne_sumvec(dim,vectorofpart,laenge,sumvec);
149 	for (i=1L;i<dim-2L;i++)
150 	   for (j=i+1L;j<dim-1L;j++)
151 	      if ( S_V_II(laenge,i) > S_V_II(laenge,j) )
152 		 m_i_i(0L,S_M_IJ(Dominanz,i,j));
153 	      else
154 		 {
155 		 m_i_i(1L,S_M_IJ(Dominanz,i,j));
156 		 for (k=0L;k<S_V_II(laenge,i);k++)
157 		    if (S_M_IJI(Dominanz,i,j) == 1L)
158 		       if (S_V_II(S_V_I(sumvec,i),k) -
159 			   S_V_II(S_V_I(sumvec,j),k) < 0 )
160 			  m_i_i(0L,S_M_IJ(Dominanz,i,j));
161 		 } /* else */
162 
163 	for (i=0L;i<dim;i++)
164 		for (j=0L;j<dim;j++)
165 				if (emptyp(S_M_IJ(Dominanz,i,j)))
166 					m_i_i(0L,S_M_IJ(Dominanz,i,j));
167 
168 	freeall(sumvec); freeall(laenge);
169 	return OK;
170 }                                        /* Ende berechne_Dominanz */
171 #endif /* PARTTRUE */
172 #endif /* MATRIXTRUE */
173 /*-----------------------------------------------------------------*/
174 /*-----------------------------------------------------------------*/
175 
176 
177 
178 
179 /*-----------------------------------------------------------------*/
180 /*-----------------------------------------------------------------*/
181 /*             siehe Programmbeschreibung unter 3.                 */
182 /*-----------------------------------------------------------------*/
183 #ifdef MATRIXTRUE
184 #ifdef PARTTRUE
berechne_sprod_kk(dim,n,vectorofpart,ytinvers,M_ks,Dominanz,sprod_kk)185 static INT berechne_sprod_kk(dim,n,vectorofpart,ytinvers,M_ks,Dominanz,
186 		     sprod_kk)
187         /* berechnet die Matrix sprod_kk, die an der Stelle (i,j) das
188            "alpha-Skalarprodukt" der monomialsym. Polynome zu den
189            Partitionen mit Indizes i und j enthaelt */
190 	INT dim;
191 	OP n,vectorofpart,ytinvers,M_ks,Dominanz,sprod_kk;
192 /* AK 210891 V1.3 */
193 {
194 	INT i,j,k;
195         INT erg = OK;
196 	OP part,alpha,l,zen,h_eins,h_zwei,zen_zwei;
197 
198 	l=callocobject();
199 	zen=callocobject();
200 	h_eins=callocobject();
201 	alpha=callocobject();
202 	h_zwei=callocobject();
203 	zen_zwei=callocobject();
204 
205 	m_il_v(dim,zen_zwei);
206 	for(i=0L;i<dim;i++)
207 	   m_i_i(0L,S_V_I(zen_zwei,i));
208 	m_i_i(2L,alpha);
209 	for(i=0L;i<dim;i++)
210 	   {
211 	    part=callocobject();
212 	    copy(S_V_I(vectorofpart,i),part);
213 	    ordcen(part,zen);
214 	    length_partition(part,l);
215 	    hoch(alpha,l,h_eins);
216 	    mult(zen,h_eins,S_V_I(zen_zwei,i));
217 	    freeall(part);
218             FREESELF(l);
219 	   }
220 	for(i=0L;i<dim;i++)
221 	   for(j=i;j<dim;j++)
222 	      {
223 	      m_i_i(0L,S_M_IJ(sprod_kk,i,j));
224 	      if(S_M_IJI(Dominanz,i,j)==1L)
225 		 for(k=0L;k<dim;k++)
226 		    {
227 		    mult(S_M_IJ(M_ks,i,k),S_M_IJ(M_ks,j,k),h_zwei);
228 		    mult(h_zwei,S_V_I(zen_zwei,k),h_zwei);
229 		    add_apply(h_zwei,S_M_IJ(sprod_kk,i,j));
230 		    }
231 	       }
232 
233 	freeall(alpha); freeall(l); freeall(zen);
234 	freeall(h_eins); freeall(h_zwei); freeall(zen_zwei);
235 	ENDR("berechne_sprod_kk:internal function");
236 }                                        /* Ende berechne_sprod_kk */
237 #endif /* PARTTRUE */
238 #endif /* MATRIXTRUE */
239 /*-----------------------------------------------------------------*/
240 /*-----------------------------------------------------------------*/
241 
242 
243 
244 
245 #ifdef MATRIXTRUE
berechne_Ainvers_und_b(dim,i,Ainvers,b,M_Pk,sprod_kk,Dominanz)246 static INT berechne_Ainvers_und_b(dim,i,Ainvers,b,M_Pk,sprod_kk,Dominanz)
247 	INT dim,i;
248 	OP Ainvers,b,M_Pk,sprod_kk,Dominanz;
249 /* AK 210891 V1.3 */
250 {
251 	INT j,k;
252 	OP h,Ainvers_alt;
253 
254 	h=callocobject();
255 	Ainvers_alt=callocobject();
256 
257 	m_ilih_m(dim-i-2L,dim-i-2L,Ainvers_alt);
258 	for (j=dim-i-3L;j>=0L;j--)
259 	   for (k=0L;k<=dim-i-3L;k++)
260 	      copy(S_M_IJ(Ainvers,j,k),S_M_IJ(Ainvers_alt,j,k));
261 
262 	m_ilih_m(dim-i-1L,dim-i-1L,Ainvers); /* Speicher fuer neues Ainvers */
263 	for (j=dim-i-3L;j>=0L;j--)
264 	   for (k=0L;k<=dim-i-3L;k++)
265 	      copy(S_M_IJ(Ainvers_alt,j,k),S_M_IJ(Ainvers,j+1L,k+1L));
266 	freeall(Ainvers_alt); /* Ainvers_alt freigeben, Werte auf neues
267 				 Ainvers kopiert */
268 	m_i_i(0L,S_M_IJ(Ainvers,0L,0L));
269 
270 	for (j=0L;j<dim;j++)
271 		for (k=0L;k<dim;k++)
272 			{
273 			if (emptyp(S_M_IJ(sprod_kk,j,k)))
274 				m_i_i(0L,S_M_IJ(sprod_kk,j,k));
275 			if (emptyp(S_M_IJ(M_Pk,j,k)))
276 				m_i_i(0L,S_M_IJ(M_Pk,j,k));
277 			}
278 
279 	for (k=i+1L;k<dim;k++)
280 	   {
281 	   mult(S_M_IJ(M_Pk,i+1L,k),S_M_IJ(sprod_kk,i+1L,k),h);
282 	   add(S_M_IJ(Ainvers,0L,0L),h,S_M_IJ(Ainvers,0L,0L));
283 	   }
284 	invers(S_M_IJ(Ainvers,0L,0L),S_M_IJ(Ainvers,0L,0L));
285 	for (j=1L;j<dim-i-1L;j++)
286 	   {
287 	   m_i_i(0L,S_M_IJ(Ainvers,j,0L));
288 	   for (k=1L;k<=j;k++)
289 	      {
290 	      mult(S_M_IJ(b,k-1L,0L),S_M_IJ(Ainvers,j,k),h);
291 	      add_apply(h,S_M_IJ(Ainvers,j,0L));
292 	      } /* for k */
293 	   mult_apply(S_M_IJ(Ainvers,0L,0L),S_M_IJ(Ainvers,j,0L));
294 	   } /* for j */
295 						  /* Ainvers ist berechnet */
296 
297 	m_ilih_m(1L,dim-i-1L,b);
298 	for (j=0L;j<dim-i-1L;j++)
299 	   {
300 	   m_i_i(0L,S_M_IJ(b,j,0L));
301 	   if(S_M_IJI(Dominanz,i,j+i+1L) == 1L)
302 	      {
303 	      for (k=j;k<dim;k++)
304 		 {
305 		 mult(S_M_IJ(M_Pk,j+i+1L,k),S_M_IJ(sprod_kk,i,k),h);
306 		 add(S_M_IJ(b,j,0L),h,S_M_IJ(b,j,0L));
307 		 } /* for k */
308 	      addinvers(S_M_IJ(b,j,0L),S_M_IJ(b,j,0L));
309 	      }
310 	   } /* for j */
311 
312 	freeall(h);
313 	return OK;
314 }                                   /* Ende berechne_Ainvers_und_b */
315 #endif /* MATRIXTRUE */
316 /*-----------------------------------------------------------------*/
317 
318 
319 
320 
321 /*-----------------------------------------------------------------*/
322 /*-----------------------------------------------------------------*/
323 /*             siehe Programmbeschreibung unter 4.                 */
324 /*-----------------------------------------------------------------*/
325 #ifdef MATRIXTRUE
berechne_M_Pk(dim,sprod_kk,vectorofpart,Dominanz,M_Pk)326 static INT berechne_M_Pk(dim,sprod_kk,vectorofpart,Dominanz,M_Pk) INT dim;
327 	OP sprod_kk,vectorofpart,Dominanz,M_Pk;
328 /* AK 210891 V1.3 */
329 {
330 	INT i,j,k,domdim,sum,index[77];
331 	OP Ainvers,b,C,Adachinv,bdach;
332 
333 	Ainvers=callocobject();
334 	b=callocobject();
335 	C=callocobject();
336 
337 	for(i=0L;i<dim;i++) m_i_i(1L,S_M_IJ(M_Pk,i,i));
338 	m_ilih_m(1L,1L,Ainvers);
339 	invers(S_M_IJ(sprod_kk,dim-1L,dim-1L),S_M_IJ(Ainvers,0L,0L));
340 	m_ilih_m(1L,1L,b);
341 	addinvers(S_M_IJ(sprod_kk,dim-2L,dim-1L),S_M_IJ(b,0L,0L));
342 	mult(Ainvers,b,C);
343 	copy(S_M_IJ(C,0L,0L),S_M_IJ(M_Pk,dim-2L,dim-1L));
344 	freeall(C);
345 
346 	for(i=dim-3L;i>=0L;i--)  /* Berechnung des zonalen Polynoms zur
347 				  Partition mit Index i */
348 	  {
349 	   C=callocobject();
350 	   Adachinv=callocobject();
351 	   bdach=callocobject();
352 
353 	   domdim=0L;
354 	   for (j=i+1L;j<dim;j++)
355 	      domdim = domdim + S_M_IJI(Dominanz,i,j);
356 				  /* domdim ist die Dimension von Adachinv */
357 	   for (j=0L;j<domdim;j++)
358 	      {
359 	      sum = 0L;
360 	      k = i;
361 	      while (sum < j+1L)
362 		 sum = sum + S_M_IJI(Dominanz,i,++k);
363 	      index[j] = k;
364 	      }
365 	   berechne_Ainvers_und_b(dim,i,Ainvers,b,M_Pk,sprod_kk,Dominanz);
366 	   m_ilih_m(domdim,domdim,Adachinv);
367 	   m_ilih_m(1L,domdim,bdach);
368 
369 	   for (k=0L;k<domdim;k++)
370 	      copy(S_M_IJ(b,index[k]-(i+1L),0L),S_M_IJ(bdach,k,0L));
371 	   for (j=0L;j<domdim;j++)
372 	      for (k=0L;k<domdim;k++)
373 		    if (k<=j)
374 			copy(S_M_IJ(Ainvers,index[j]-(i+1L),index[k]-(i+1L)),
375 			     S_M_IJ(Adachinv,j,k));
376 	   for (j=0L;j<domdim;j++)
377 	      for(k=j+1L;k<domdim;k++)
378 		 m_i_i(0L,S_M_IJ(Adachinv,j,k));
379 	   mult_udrmatrix_matrix(Adachinv,bdach,C);
380 	    /* Multiplizieren der unteren Dreieckmatrix Adachinv mit bdach */
381 
382 	   for (j=0L;j<domdim;j++)
383 	      copy(S_M_IJ(C,j,0L),S_M_IJ(M_Pk,i,index[j]));
384 	   for (j=i+1L;j<index[0];j++)
385 	      m_i_i(0L,S_M_IJ(M_Pk,i,j));
386 	   for (j=0L;j<domdim-1L;j++)
387 	      for (k=index[j]+1L;k<index[j+1];k++)
388 		 m_i_i(0L,S_M_IJ(M_Pk,i,k));
389 
390 	   freeall(Adachinv);
391 	   freeall(bdach);
392 	   freeall(C);
393 	  }
394 	freeall(Ainvers);
395 	freeall(b);
396 	return OK;
397 }                                            /* Ende berechne_M_Pk */
398 #endif /* MATRIXTRUE */
399 /*-----------------------------------------------------------------*/
400 /*-----------------------------------------------------------------*/
401 
402 
403 
404 
405 /*-----------------------------------------------------------------*/
406 /*-----------------------------------------------------------------*/
407 /*             siehe Programmbeschreibung unter 5.                 */
408 /*-----------------------------------------------------------------*/
409 #ifdef KOSTKATRUE
berechne_M_ke(dim,n,K,vectorofpart,M_ke)410 static INT berechne_M_ke(dim,n,K,vectorofpart,M_ke)
411      /* M_ke ist die Uebergangsmatrix zwischen den monomialsymm.
412         Polynomen und den Produkten der elementarsymm. Funktionen. */
413 	INT dim; OP n,K,vectorofpart,M_ke;
414 /* AK 210891 V1.3 */
415 {
416 	OP J,Kinvers,Ktinv;
417 
418 	Kinvers=callocobject();
419 	Ktinv=callocobject();
420 	J=callocobject();
421 
422 	invers_matrix(K,Kinvers);
423 	transpose_matrix(Kinvers,Ktinv);
424 	make_n_transpositionmatrix(n,J);
425 	mult(Ktinv,J,M_ke);
426 	mult(M_ke,Kinvers,M_ke);
427 
428 	freeall(J); freeall(Kinvers); freeall(Ktinv);
429 	return OK;
430 }/* Ende berechne_M_ke */
431 #endif /* KOSTKATRUE */
432 /*-----------------------------------------------------------------*/
433 /*-----------------------------------------------------------------*/
434 
435 
436 
437 
438 /*-----------------------------------------------------------------*/
439 /*-----------------------------------------------------------------*/
440 /*             siehe Programmbeschreibung unter 6.                 */
441 /*-----------------------------------------------------------------*/
442 #ifdef MATRIXTRUE
Anpassen(dim,n,ytstern,M_ke,M_Pk)443 static INT Anpassen(dim,n,ytstern,M_ke,M_Pk)
444            /* M_Pk wird anders normiert (siehe James, z.B. 1964), man
445               erhaelt M_Zk. Daraus wird M_Zs, die Uebergangsmatrix
446               zwischen den zonalen Polynomen Z und den Produkten der
447               Potenzsummen berechnet. Diese M_Zs sind tabelliert.
448               Analoges gilt fuer die M_Ze. */
449 	INT dim;
450 	OP n,ytstern,M_ke,M_Pk;
451 /* AK 210891 V1.3 */
452 {
453 	INT i,j;
454 	OP h,diag,M_Zk,M_Zs,M_Ze;
455 	h=callocobject();
456 	diag=callocobject();
457 	M_Zk=callocobject();
458 	M_Zs=callocobject();
459 	M_Ze=callocobject();
460 
461 	m_ilih_m(dim,dim,diag);
462 	fakul(n,h);
463 	for(i=0L;i<dim;i++)
464 	  {
465 	   for(j=0L;j<i;j++)
466 	     {
467 	      m_i_i(0L,S_M_IJ(M_Pk,i,j));
468 	      m_i_i(0L,S_M_IJ(diag,i,j));
469 	     }
470 	   for(j=i;j<dim;j++)
471 	      m_i_i(0L,S_M_IJ(diag,i,j));
472 	  }
473 	for(i=0L;i<dim;i++)
474 	   div(h,S_M_IJ(M_Pk,i,dim-1L),S_M_IJ(diag,i,i));
475 	m_ilih_m(dim,dim,M_Zk);
476 	mult(diag,M_Pk,M_Zk);
477 
478 	printf("Uebergangsmatrizen, um die zonalen Polynome Z in \n");
479 	printf("verschiedenen Basen darzustellen:\n");
480 	printf("\n");
481 	printf("Basis : monomialsymm. Funktionen\n");
482 	println(M_Zk);
483 	printf("\n");
484 	printf("Basis : Potenzsummen\n");
485 	m_ilih_m(dim,dim,M_Zs);
486 	mult(M_Zk,ytstern,M_Zs);
487 	println(M_Zs);
488 	printf("\n");
489 	printf("Basis : elementarsymm. Funktionen\n");
490 	m_ilih_m(dim,dim,M_Ze);
491 	mult(M_Zk,M_ke,M_Ze);
492 	println(M_Ze);
493 	printf("\n");
494 	printf("\n");
495 
496 
497 	freeall(h); freeall(diag); freeall(M_Zk); freeall(M_Zs);
498 	freeall(M_Ze);
499 	return OK;
500 }                                                 /* Ende Anpassen */
501 
502 
503 
mult_udrmatrix_matrix(a,b,c)504 static INT mult_udrmatrix_matrix(a,b,c) OP a,b,c;
505 /* AK 280388 */
506 /* AK 060789 V1.0 */ /* AK 210891 V1.3 */
507 {
508 	INT i,j,k;
509         INT erg = OK;
510 	OP z; /* zwischen ergebnis bei matrix-multiplikation */
511 
512 	if (neq(s_m_l(a),s_m_h(b)))
513 	{
514 		error("mult_matrix_matrix:can not be multiplied");
515 		return(ERROR);
516 	};
517 
518 	m_ilih_nm(S_M_LI(b),S_M_HI(a),c);
519 	z=callocobject(); /* zwischensumme*/
520 	for (i=0L;i<S_M_HI(a);i++)	/* ueber zeilen der linken Matrix */
521 		for (j=0L;j<S_M_LI(b);j++) /* ueber spalten der rechten Matrix */
522 			for (k=0L;k<i+1L;k++)
523 			{
524 				MULT(S_M_IJ(a,i,k),S_M_IJ(b,k,j),z);
525 				add_apply(z,S_M_IJ(c,i,j));
526                                 FREESELF(z);
527 			};
528 	FREEALL(z);
529 	ENDR("mult_udrmatrix_matrix:internal function");
530 }
531 #endif
532 
533 
534