1 /* multiplication homsym \times monomial = monomial */
2 #include "def.h"
3 #include "macro.h"
4
5
6
7 INT mhm_integer_partition_();
8 INT mhm_integer_partition_hashtable();
9 INT mhm_integer_hashtable_hashtable();
10 INT thm_integer__faktor();
11
12 static INT SYMMETRICA_mhm_co_ip();
13 static INT hm_coeff();
14
mhm_null__(b,c,f)15 INT mhm_null__(b,c,f) OP b,c,f;
16 {
17 INT erg = OK;
18 CTTTO(HASHTABLE,PARTITION,MONOMIAL,"mhm_null__(1)",b);
19 CTO(HASHTABLE,"mhm_null__(2)",c);
20 if (S_O_K(b) == PARTITION) {
21 _NULL_PARTITION_(b,c,f);
22 }
23 else /* monomial or hashtable */
24 {
25 OP z;
26 FORALL(z,b,{
27 OP d;
28 d = CALLOCOBJECT();
29 erg += b_sk_mo(CALLOCOBJECT(),CALLOCOBJECT(),d);
30 erg += copy_partition(S_MO_S(z),S_MO_S(d));
31 COPY(S_MO_K(z),S_MO_K(d));
32 if (not EINSP(f))
33 {
34 MULT_APPLY(f,S_MO_K(d));
35 }
36 INSERT_HASHTABLE(d,c,add_koeff,eq_monomsymfunc,hash_monompartition);
37 });
38 }
39 ENDR("mhm_null__");
40 }
41
mhm_integer__(a,b,c,f)42 INT mhm_integer__(a,b,c,f) OP a,b,c; OP f;
43 /* AK 311001 */
44 {
45 INT erg = OK;
46
47 CTO(INTEGER,"mhm_integer__(1)",a);
48 CTTTO(HASHTABLE,PARTITION,MONOMIAL,"mhm_integer__(2)",b);
49 CTO(HASHTABLE,"mhm_integer__(3)",c);
50 SYMCHECK((S_I_I(a) < 0),"mhm_integer__:parameter < 0");
51 if (S_I_I(a) == 0) {
52 erg += mhm_null__(b,c,f);
53 goto ende;
54 }
55 else if (S_O_K(b) == PARTITION) {
56 erg += mhm_integer_partition_hashtable(a,b,c,f);
57 goto ende;
58 }
59 else
60 {
61 erg += mhm_integer_hashtable_hashtable(a,b,c,f);
62 goto ende;
63 }
64 ende:
65 ENDR("mhm_integer__");
66 }
67
mhm_partition__(a,b,c,f)68 INT mhm_partition__(a,b,c,f) OP a,b,c; OP f;
69 /* AK 311001 */
70 {
71 INT erg = OK;
72 CTO(PARTITION,"mhm_partition__(1)",a);
73 CTTTO(HASHTABLE,PARTITION,MONOMIAL,"mhm_partition__(2)",b);
74 CTO(HASHTABLE,"mhm_partition__(3)",c);
75
76 if (S_PA_LI(a) == 0) {
77 erg += mhm_null__(b,c,f);
78 goto ende;
79 }
80 else if (S_PA_LI(a) == 1) {
81 erg += mhm_integer__(S_PA_I(a,0),b,c,f);
82 goto ende;
83 }
84 else { /* partition of length >= 1 */
85 INT i; OP d,e;
86
87 d=CALLOCOBJECT();
88 e=CALLOCOBJECT();
89
90 erg += init_hashtable(e);
91 erg += thm_integer__faktor(S_PA_I(a,0),e,f);
92 for (i=1;i<S_PA_LI(a);i++)
93 {
94 FREESELF(d);
95 erg += init_hashtable(d);
96 SWAP(d,e);
97 erg += mhm_integer__(S_PA_I(a,i),d,e,cons_eins);
98 }
99 FREEALL(d);
100
101 mult_monomial_monomial(e,b,c);
102 FREEALL(e);
103 }
104
105 ende:
106 ENDR("mhm_partition__");
107 }
108
mhm_homsym__(a,b,c,f)109 INT mhm_homsym__(a,b,c,f) OP a,b,c,f;
110 /* AK 061101 */
111 /* c += h_a \times s_b \times f */
112 {
113 INT erg = OK;
114 OP z,ff;
115 CTO(HOMSYM,"mhm_homsym__(1)",a);
116 CTTTO(HASHTABLE,PARTITION,MONOMIAL,"mhm_homsym__(2)",b);
117 CTO(HASHTABLE,"mhm_homsym__(3)",c);
118
119 if (S_L_S(a) == NULL) { /* empty list */
120 goto eee;
121 }
122 if (S_S_N(a) == NULL) {
123 if (EINSP(f))
124 erg += mhm_partition__(S_S_S(a),b,c,S_S_K(a));
125 else {
126 ff = CALLOCOBJECT();
127 MULT(f,S_S_K(a),ff);
128 erg += mhm_partition__(S_S_S(a),b,c,ff);
129 FREEALL(ff);
130 }
131 goto eee;
132 }
133 else if (S_S_SLI(a) == 0) {
134 if (EINSP(f))
135 erg += mhm_partition__(S_S_S(a),b,c,S_S_K(a));
136 else {
137 ff = CALLOCOBJECT();
138 MULT(f,S_S_K(a),ff);
139 erg += mhm_partition__(S_S_S(a),b,c,ff);
140 FREEALL(ff);
141 }
142 erg += mhm_homsym__(S_S_N(a),b,c,f);
143 goto eee;
144 }
145 else {
146 OP zv,hi;
147 INT i,j;
148 i = S_S_SII(a,0);
149 z = a;zv = NULL;
150 again:
151 if (S_S_SII(z,0) == i)
152 {
153 for (j=0;j<S_S_SLI(z)-1;j++)
154 M_I_I(S_S_SII(z,j+1),S_S_SI(z,j));
155 M_I_I(S_S_SLI(z)-1,S_S_SL(z));
156 zv = z; z = S_S_N(z);
157 if (z != NULL) goto again;
158 }
159 /* jetzt ist z der erste teil mit dem kleinsten teil partition >i */
160 /* zv ist der letzte teil mit kleinsten teil = i */
161 /* berechne : h_i * b * h_a + b* h_z */
162 C_S_N(zv,NULL);
163
164 ff = CALLOCOBJECT();
165 init_hashtable(ff);
166 hi = CALLOCOBJECT();
167 M_I_I(i,hi);
168 erg += mhm_integer__(hi,b,ff,f);
169 erg += mhm_homsym__(a,ff,c,cons_eins);
170
171 if (z != NULL)
172 erg += mhm_homsym__(z,b,c,f);
173 FREEALL(hi);
174 FREEALL(ff);
175 /* a wieder richtig zusammen bauen */
176
177 zv = a;
178 aa:
179 for (j=S_S_SLI(zv);j>0;j--)
180 M_I_I(S_S_SII(zv,j-1),S_S_SI(zv,j));
181 M_I_I(i,S_S_SI(zv,0));
182 M_I_I(S_S_SLI(zv)+1,S_S_SL(zv));
183 if (S_S_N(zv) != NULL) { zv = S_S_N(zv); goto aa; }
184
185 C_S_N(zv,z);
186
187 }
188
189
190 eee:
191 ENDR("mhm_homsym__");
192 }
193
mhm_hashtable__(a,b,c,f)194 INT mhm_hashtable__(a,b,c,f) OP a,b,c,f;
195 /* AK 061101 */
196 /* c += h_a \times s_b \times f */
197 {
198 INT erg = OK;
199 CTO(HASHTABLE,"mhm_hashtable__(1)",a);
200 CTTTO(HASHTABLE,PARTITION,MONOMIAL,"mhm_hashtable__(2)",b);
201 CTO(HASHTABLE,"mhm_hashtable__(3)",c);
202 M_FORALL_MONOMIALS_IN_A(a,b,c,f,mhm_partition__);
203 ENDR("mhm_homsym__");
204 }
205
mhm_integer_partition_hashtable(a,b,c,f)206 INT mhm_integer_partition_hashtable(a,b,c,f) OP a,b,c,f;
207 /* AK 061101 */
208 {
209 INT erg = OK;
210
211 CTO(INTEGER,"mhm_integer_partition_(1)",a);
212 CTO(PARTITION,"mhm_integer_partition_(2)",b);
213 CTO(HASHTABLE,"mhm_integer_partition_(3)",c);
214 SYMCHECK((S_I_I(a) < 0),"mhm_integer_partition_hashtable:integer < 0");
215
216 if (S_I_I(a) == 0) {
217 erg += mhm_null__(b,c,f);
218 goto ende;
219 }
220 else if (S_PA_LI(b) == 0)
221 {
222 erg += thm_integer__faktor(a,c,f); /* generates sum over all partitions */
223 goto ende;
224 }
225 else
226 {
227 erg += SYMMETRICA_mhm_co_ip(a,b,c,f);
228 goto ende;
229 }
230
231 ende:
232 ENDR("mhm_integer_partition_hashtable");
233 }
234
mhm_integer_hashtable_hashtable(a,b,c,f)235 INT mhm_integer_hashtable_hashtable(a,b,c,f) OP a,b,c,f;
236 /* AK 061101 */
237 {
238 INT erg = OK;
239 CTO(INTEGER,"mhs_integer_hashtable_(1)",a);
240 CTTO(HASHTABLE,MONOMIAL,"mhs_integer_hashtable_(2)",b);
241 CTO(HASHTABLE,"integer_hashtable_(3)",c);
242 M_FORALL_MONOMIALS_IN_B(a,b,c,f,mhm_integer_partition_hashtable);
243 ENDR("mhm_integer_hashtable_hashtable");
244 }
245
246
247
mult_homsym_monomial(a,b,c)248 INT mult_homsym_monomial(a,b,c) OP a,b,c;
249 /* AK 111001
250 */
251 {
252 INT erg = OK;
253 INT t=0; /* is 1 if transfer HASHTABLE->MONOMIAL necessary */
254 CTTTTO(HASHTABLE,INTEGER,PARTITION,HOMSYM,"mult_homsym_monomial(1)",a);
255 CTTTO(HASHTABLE,PARTITION,MONOMIAL,"mult_homsym_monomial(2)",b);
256 CTTTO(EMPTY,HASHTABLE,MONOMIAL,"mult_homsym_monomial(3)",c);
257
258 if (S_O_K(c) == EMPTY) {
259 t=1;
260 init_hashtable(c);
261 }
262 else if (S_O_K(c) == MONOMIAL) {
263 t=1;
264 t_MONOMIAL_HASHTABLE(c,c);
265 }
266
267 if (S_O_K(a) == INTEGER)
268 {
269 erg += mhm_integer__(a,b,c,cons_eins);
270 goto ende;
271 }
272 else if (S_O_K(a) == PARTITION)
273 {
274 erg += mhm_partition__(a,b,c,cons_eins);
275 goto ende;
276 }
277 else if (S_O_K(a) == HOMSYM)
278 {
279 erg += mhm_homsym__(a,b,c,cons_eins);
280 goto ende;
281 }
282 else /* if (S_O_K(a) == HASHTABLE) */
283 {
284 erg += mhm_hashtable__(a,b,c,cons_eins);
285 goto ende;
286 }
287
288 ende:
289 if (t==1) t_HASHTABLE_MONOMIAL(c,c);
290 ENDR("mult_homsym_monomial");
291 }
292
next_part_EXPONENT_apply_limit(p,be,bp)293 static INT next_part_EXPONENT_apply_limit(p,be,bp) OP p,be; INT bp;
294 /* next with limit be */
295 /* bp maximaler eintrag != 0 */
296 {
297 INT j,i,w,t,h;
298 INT erg = OK;
299 OP z,zz;
300
301 i = S_PA_LI(p);t=0;
302 for (j=bp, z = S_PA_I(p,bp);j>0;j--,z--)
303 {
304 w = S_I_I(z);
305 if (w > 0) /* schauen ob frei */
306 {
307 t+=w;
308 if (j >= S_PA_LI(be))
309 {
310 i = j;
311 }
312 else if (S_PA_II(be,j) == 0)
313 {
314 i = j;
315 }
316 else /* S_PA_II(be,j) > 0 */ {
317 t -= S_PA_II(be,j);
318 if (t>0)
319 i = j;
320 }
321 }
322 else {
323 if (j < S_PA_LI(be))
324 {
325 if (S_PA_II(be,j) > 0) t-=S_PA_II(be,j);
326 }
327 }
328 }
329
330 if (i == S_PA_LI(p)) return LASTPARTITION;
331 /* an der stelle i kann decrementiert werden */
332 w=0;
333 for(j=1,z=S_V_S(S_PA_S(p));j<=i;j++,z++)
334 {
335 h = S_I_I(z);
336 C_I_I(z,0);
337 while (h--) w+=j;
338 }
339
340 w += (i+1);
341 DEC_INTEGER(S_PA_I(p,i));
342 /* nun das gewicht w in den kleineren teilen unter bringen,
343 aber das limit beruecksichtigen */
344
345 /* zuerst feststellen wie viel vom limit bereits abgearbeitet ist */
346 for (t=0,j=bp;j>=i;j--)
347 t+= S_PA_II(p,j);
348 /* t teile in p > i */
349 /* nun feststellen wieviel teile davon das limit ueberdecken */
350 for (j=S_PA_LI(be)-1,z = S_PA_I(p,j),
351 zz = S_PA_I(be,j);j>=0;j--,z--,zz--)
352 {
353 h = S_I_I(zz);
354 if (t >= h)
355 t -= h;
356 else if (t > 0)
357 {
358 C_I_I(z,h-t);
359 t = 0;
360 w -= S_I_I(z)*(j+1);
361 }
362 else /* t == 0 */
363 {
364 C_I_I(z,h);
365 w -= h*(j+1);
366 }
367 }
368 i--;
369 t = i;
370 while (w > 0) {
371 if (i==0) /* kein limit mehr, bei t einfuegen */
372 {
373 M_I_I( S_PA_II(p,t) + (w/ (t+1)) ,S_PA_I(p,t));
374 if ( (w % (t+1) ) > 0 )
375 INC_INTEGER(S_PA_I(p, w % (t+1) -1));
376 return OK;
377 }
378
379 /* nach links kleineres teil aus dem limit suchen und vergroessern */
380 for (j=i-1, z = S_PA_I(p,j);j>=0;j--,z--)
381 {
382 nochmal:
383 if ((i-j) > w) continue;
384 if (S_I_I(z) > 0) {
385 DEC_INTEGER(z);
386 INC_INTEGER(S_PA_I(p,i));
387 w += j;
388 w -= i;
389 if (w == 0) return OK;
390 goto nochmal;
391 }
392 }
393 i--;
394 }
395 SYMCHECK(1,"next_part_EXPONENT_apply_limit:should never be here");
396 ENDR("next_part_EXPONENT_apply_limit");
397 }
398
399 static OP hm_coeff_lo,hm_coeff_oben,hm_coeff_unten;
SYMMETRICA_mhm_co_ip(a,b,c,faktor)400 static INT SYMMETRICA_mhm_co_ip(a,b,c,faktor) OP a,b,c; OP faktor;
401 /* addiert das ergebnis von h_a \times m_b zu c */
402 /* dabei gibt es den coeffizienten faktor */
403 {
404 INT erg = OK;
405 INT bp;
406 INT i,w,j;
407 OP p;
408 OP be,m,z;
409 CTO(INTEGER,"SYMMETRICA_mhm_co_ip(1)",a);
410 CTO(PARTITION,"SYMMETRICA_mhm_co_ip(2)",b);
411 CTO(HASHTABLE,"SYMMETRICA_mhm_co_ip(3)",c);
412
413 hm_coeff_lo = CALLOCOBJECT();
414 hm_coeff_oben=CALLOCOBJECT(); C_O_K(hm_coeff_oben,INTEGER);
415 hm_coeff_unten=CALLOCOBJECT(); C_O_K(hm_coeff_unten,INTEGER);
416
417 for (i=0,w=S_I_I(a),z=S_V_S(S_PA_S(b));i<S_PA_LI(b);i++,z++) w+= S_I_I(z);
418 /* w is the weight of the result */
419
420
421 p=CALLOCOBJECT();
422 be=CALLOCOBJECT();
423
424
425 erg += t_VECTOR_EXPONENT(b,be);
426
427 erg += copy_partition(be,p);
428 /* increase the partition, filling with zero */
429 i = S_PA_LI(p);
430
431
432
433 inc_vector_co(S_PA_S(p),S_I_I(a));
434 for (;i<S_PA_LI(p);i++) M_I_I(0,S_PA_I(p,i));
435
436
437
438 for (i=S_PA_LI(p)-1;i>=0;i--)
439 if (S_PA_II(p,i) > 0 ) break;
440 DEC_INTEGER(S_PA_I(p,i));
441 M_I_I(1,S_PA_I(p,i+S_I_I(a)));
442
443 m = CALLOCOBJECT();
444 erg += b_sk_mo(CALLOCOBJECT(),CALLOCOBJECT(),m);
445 erg += b_ks_pa(VECTOR,CALLOCOBJECT(),S_MO_S(m));
446 erg += m_il_nv(w,S_PA_S(S_MO_S(m)));
447 C_O_K(S_PA_S(S_MO_S(m)), INTEGERVECTOR);
448
449 do {
450 FREESELF(S_MO_K(m));
451 erg += hm_coeff(be,p,S_MO_K(m));
452 {
453 INT i,j=0,ba=0;
454 OP l;
455 for (l=S_V_S(S_PA_S(p)),i=0; i<S_PA_LI(p); i++,l++)
456 if (S_I_I(l)>0)
457 {
458 j += S_I_I(l);
459 ba=i;
460 }
461
462 /* ba is the last non zero entry in p */
463 C_I_I(S_PA_L(S_MO_S(m)),j);
464
465 for (l=S_V_S(S_PA_S(S_MO_S(m))),i=0;i<=ba;i++)
466 if (S_PA_II(p,i)>0)
467 for (j=(INT)0;j<S_PA_II(p,i);j++)
468 {
469 C_I_I(l,i+1);
470 l++;
471 }
472 }
473
474 bp = S_PA_II(S_MO_S(m), S_PA_LI(S_MO_S(m))-1) -1;
475 HASH_INTEGERVECTOR(S_PA_S(S_MO_S(m)),j);
476 C_PA_HASH(S_MO_S(m),j);
477 if (not EINSP(faktor))
478 MULT_APPLY(faktor,S_MO_K(m));
479
480 add_apply_hashtable(m,c,add_koeff,eq_monomsymfunc,hash_monompartition);
481
482 } while(next_part_EXPONENT_apply_limit(p,be,bp) != LASTPARTITION);
483
484 FREEALL(m);
485 FREEALL(p);
486 FREEALL(be);
487 FREEALL(hm_coeff_lo);
488 FREEALL(hm_coeff_oben);
489 FREEALL(hm_coeff_unten);
490 ENDR("internal to mhm(2)");
491 }
492
493
494
hm_coeff(ae,be,c)495 static INT hm_coeff(ae,be,c) OP ae,be,c;
496 /* ae ist kleine partition (exponent), be ist grosse partition (exponent)*/
497 /* c wird das ergebnis, der coeff von be beim produkt h_n * ae */
498 {
499 INT i,j,k,t=0,w;
500 INT erg = OK;
501 OP z;
502 CTO(EMPTY,"hm_coeff(3)",c);
503 M_I_I(1,c);
504
505 for (i=S_PA_LI(ae)-1,j=S_PA_LI(be)-1,k=0,
506 z = S_PA_I(be,j);
507 i>=0;i--)
508 {
509 w = S_PA_II(ae,i);
510 if (w > 0) {
511 while(j>=i) {
512 k+= S_I_I(z);
513 j--;z--;
514 }
515 if (k < w)
516 {
517 if (S_O_K(c) == INTEGER) M_I_I(0,c);
518 else m_i_i(0,c);
519 goto endr_ende;
520 }
521
522 C_I_I(hm_coeff_oben,k);
523 C_I_I(hm_coeff_unten,w);
524 if (t==0)
525 {
526 BINOM_POSINTEGER_POSINTEGER(hm_coeff_oben,hm_coeff_unten,c);
527 t=1;
528 }
529 else {
530 BINOM_POSINTEGER_POSINTEGER(hm_coeff_oben,hm_coeff_unten,hm_coeff_lo);
531 MULT_APPLY(hm_coeff_lo,c);
532 FREESELF(hm_coeff_lo);
533 }
534
535 k -= w;
536 }
537 }
538
539
540 ENDR("internal(2) to mult_homsym_monomial");
541 }
542
543