1 /* Littlewood-Richardson Calculator
2 * Copyright (C) 1999- Anders S. Buch (asbuch at math rutgers edu)
3 * See the file LICENSE for license information.
4 */
5
6 #include "vector.h"
7 #include "hashtab.h"
8 #include "claim.h"
9
10 #include "lincomb.h"
11 #include "schublib.h"
12
trans(vector * w,int vars,hashtab * res)13 hashtab *trans(vector *w, int vars, hashtab *res)
14 {
15 hashtab *tmp;
16 hash_itr itr;
17 vector *v;
18 int n, nw, r, s, wr, ws, i, vi, vr, last;
19
20 if (res == NULL)
21 res = hash_new((cmp_t) v_cmp, (hash_t) v_hash);
22 else
23 hash_reset(res);
24
25 nw = v_length(w);
26 n = nw;
27 while (n > 1 && v_elem(w, n-1) == n)
28 n--;
29 v_length(w) = n;
30
31 r = n-1;
32 while (r > 0 && v_elem(w, r-1) < v_elem(w, r))
33 r--;
34 if (r == 0)
35 {
36 vector *xx = v_new_zero(vars ? vars : 1);
37 hash_insertint(res, xx, 1);
38 v_length(w) = nw;
39 return res;
40 }
41 if (vars < r)
42 vars = r;
43
44 s = r + 1;
45 while (s < n && v_elem(w, r-1) > v_elem(w, s))
46 s++;
47
48 wr = v_elem(w, r-1);
49 ws = v_elem(w, s-1);
50
51 v = w;
52 v_elem(v, s-1) = wr;
53 v_elem(v, r-1) = ws;
54
55 tmp = trans(v, vars, NULL);
56 for (hash_first(tmp, itr); hash_good(itr); hash_next(itr))
57 {
58 vector *xx = hash_key(itr);
59 v_elem(xx, r-1)++;
60 hash_insertint(res, xx, hash_intvalue(itr));
61 }
62
63 last = 0;
64 vr = v_elem(v, r-1);
65 for (i = r-1; i >= 1; i--)
66 {
67 vi = v_elem(v, i-1);
68 if (last < vi && vi < vr)
69 {
70 last = vi;
71 v_elem(v, i-1) = vr;
72 v_elem(v, r-1) = vi;
73 trans(v, vars, tmp);
74 lincomb_add_multiple(res, 1, tmp, (freekey_t) v_free1, NULL);
75 v_elem(v, i-1) = vi;
76 }
77 }
78
79 v_length(w) = nw;
80 v_elem(w, s-1) = ws;
81 v_elem(w, r-1) = wr;
82 hash_free(tmp);
83
84 return res;
85 }
86
87 #if 0
88 hashtab *monk(int i, hashtab *slc, int rank)
89 {
90 hashtab *res = hash_new((cmp_t) v_cmp, (hash_t) v_hash);
91 hash_itr itr;
92
93 for (hash_first(slc, itr); hash_good(itr); hash_next(itr))
94 {
95 vector *u, *w = hash_key(itr);
96 int j, t;
97 int c = hash_intvalue(itr);
98 int n = v_length(w);
99 int wi = (i <= n) ? v_elem(w, i-1) : i;
100
101 if (i <= n+1)
102 {
103 int last, ulen = (i > n) ? i : n;
104 last = 0;
105 for (j = i-1; j >= 1; j--)
106 if (last < v_elem(w, j-1) && v_elem(w, j-1) < wi)
107 {
108 last = v_elem(w, j-1);
109 u = v_new(ulen);
110 for (t = 0; t < ulen; t++)
111 v_elem(u, t) = v_elem(w, t);
112 v_elem(u, j-1) = wi;
113 v_elem(u, i-1) = last;
114
115 if (lincomb_add_element(res, -c, u, (freekey_t) v_free1) == 0)
116 v_free(u);
117 }
118 }
119 else
120 {
121 u = v_new(i);
122 for (t = 0; t < n; t++)
123 v_elem(u, t) = v_elem(w, t);
124 for (t = n; t < i-2; t++)
125 v_elem(u, t) = t+1;
126 v_elem(u, i-2) = i;
127 v_elem(u, i-1) = i-1;
128
129 if (lincomb_add_element(res, -c, u, (freekey_t) v_free1) == 0)
130 v_free(u);
131 }
132
133 if (i >= n+1)
134 {
135 u = v_new(i + 1);
136 for (t = 0; t < n; t++)
137 v_elem(u, t) = v_elem(w, t);
138 for (t = n; t < i; t++)
139 v_elem(u, t) = t+1;
140 v_elem(u, i-1) = i+1;
141 v_elem(u, i) = i;
142
143 if (lincomb_add_element(res, c, u, (freekey_t) v_free1) == 0)
144 v_free(u);
145 }
146 else
147 {
148 int last;
149 last = (((unsigned) -1) >> 1);
150 for (j = i+1; j <= n; j++)
151 if (wi < v_elem(w, j-1) && v_elem(w, j-1) < last)
152 {
153 last = v_elem(w, j-1);
154 u = v_new(n);
155 for (t = 0; t < n; t++)
156 v_elem(u, t) = v_elem(w, t);
157 v_elem(u, i-1) = last;
158 v_elem(u, j-1) = wi;
159
160 if (lincomb_add_element(res, c, u, (freekey_t) v_free1) == 0)
161 v_free(u);
162 }
163 if (last > n && n < rank)
164 {
165 u = v_new(n+1);
166 for (t = 0; t < n; t++)
167 v_elem(u, t) = v_elem(w, t);
168 v_elem(u, i-1) = n+1;
169 v_elem(u, n) = wi;
170
171 if (lincomb_add_element(res, c, u, (freekey_t) v_free1) == 0)
172 v_free(u);
173 }
174 }
175 }
176
177 return res;
178 }
179 #endif
180
181 void _monk_add(int i, hashtab *slc, int rank, hashtab *res);
182
monk(int i,hashtab * slc,int rank)183 hashtab *monk(int i, hashtab *slc, int rank)
184 {
185 hashtab *res = hash_new((cmp_t) v_cmp, (hash_t) v_hash);
186 _monk_add(i, slc, rank, res);
187 return res;
188 }
189
_monk_add(int i,hashtab * slc,int rank,hashtab * res)190 void _monk_add(int i, hashtab *slc, int rank, hashtab *res)
191 {
192 hash_itr itr;
193
194 for (hash_first(slc, itr); hash_good(itr); hash_next(itr))
195 {
196 vector *u, *w = hash_key(itr);
197 int j, t;
198 int c = hash_intvalue(itr);
199 int n = v_length(w);
200 int wi = (i <= n) ? v_elem(w, i-1) : i;
201
202 if (i <= n+1)
203 {
204 int last, ulen = (i > n) ? i : n;
205 last = 0;
206 for (j = i-1; j >= 1; j--)
207 if (last < v_elem(w, j-1) && v_elem(w, j-1) < wi)
208 {
209 last = v_elem(w, j-1);
210 u = v_new(ulen);
211 for (t = 0; t < ulen; t++)
212 v_elem(u, t) = v_elem(w, t);
213 v_elem(u, j-1) = wi;
214 v_elem(u, i-1) = last;
215
216 if (lincomb_add_element(res, -c, u, (freekey_t) v_free1) == 0)
217 v_free(u);
218 }
219 }
220 else
221 {
222 u = v_new(i);
223 for (t = 0; t < n; t++)
224 v_elem(u, t) = v_elem(w, t);
225 for (t = n; t < i-2; t++)
226 v_elem(u, t) = t+1;
227 v_elem(u, i-2) = i;
228 v_elem(u, i-1) = i-1;
229
230 if (lincomb_add_element(res, -c, u, (freekey_t) v_free1) == 0)
231 v_free(u);
232 }
233
234 if (i >= n+1)
235 {
236 u = v_new(i + 1);
237 for (t = 0; t < n; t++)
238 v_elem(u, t) = v_elem(w, t);
239 for (t = n; t < i; t++)
240 v_elem(u, t) = t+1;
241 v_elem(u, i-1) = i+1;
242 v_elem(u, i) = i;
243
244 if (lincomb_add_element(res, c, u, (freekey_t) v_free1) == 0)
245 v_free(u);
246 }
247 else
248 {
249 int last;
250 last = (((unsigned) -1) >> 1);
251 for (j = i+1; j <= n; j++)
252 if (wi < v_elem(w, j-1) && v_elem(w, j-1) < last)
253 {
254 last = v_elem(w, j-1);
255 u = v_new(n);
256 for (t = 0; t < n; t++)
257 v_elem(u, t) = v_elem(w, t);
258 v_elem(u, i-1) = last;
259 v_elem(u, j-1) = wi;
260
261 if (lincomb_add_element(res, c, u, (freekey_t) v_free1) == 0)
262 v_free(u);
263 }
264 if (last > n && n < rank)
265 {
266 u = v_new(n+1);
267 for (t = 0; t < n; t++)
268 v_elem(u, t) = v_elem(w, t);
269 v_elem(u, i-1) = n+1;
270 v_elem(u, n) = wi;
271
272 if (lincomb_add_element(res, c, u, (freekey_t) v_free1) == 0)
273 v_free(u);
274 }
275 }
276 }
277 }
278
279
280 #if 0
281 hashtab *mult_poly_schubert(hashtab *poly, vector *perm, int rank)
282 {
283 int maxvar, i, c;
284 list *tmplist;
285 hashtab *res, *res1, *poly1;
286 hash_itr itr;
287
288 if (hash_card(poly) == 0)
289 return poly;
290
291 maxvar = 0;
292 for (hash_first(poly, itr); hash_good(itr); hash_next(itr))
293 {
294 vector *xx = hash_key(itr);
295 c = hash_intvalue(itr);
296 i = v_length(xx);
297 if (c != 0 && i > maxvar)
298 {
299 while (i > maxvar && v_elem(xx, i-1) == 0)
300 i--;
301 if (i > maxvar)
302 maxvar = i;
303 }
304 }
305
306 if (maxvar == 0)
307 {
308 hash_first(poly, itr);
309 c = hash_intvalue(itr);
310 v_free(hash_key(itr));
311 hash_reset(poly); /* LOOK OUT FOR MEMORY LEAK!! */
312 hash_insertint(poly, v_new_copy(perm), c);
313 return poly;
314 }
315
316 tmplist = l_new();
317 for (hash_first(poly, itr); hash_good(itr); hash_next(itr))
318 {
319 vector *xx = hash_key(itr);
320 if (v_length(xx) >= maxvar && v_elem(xx, maxvar-1) != 0)
321 {
322 l_append(tmplist, xx);
323 l_appendint(tmplist, hash_intvalue(itr));
324 }
325 }
326 poly1 = hash_new((cmp_t) v_cmp, (hash_t) v_hash);
327 for (i = 0; i < l_length(tmplist); i += 2)
328 {
329 vector *xx = l_elem(tmplist, i);
330 c = l_intelem(tmplist, i+1);
331 hash_remove(poly, xx);
332 v_elem(xx, maxvar-1)--;
333 while (v_length(xx) > 0 && v_elem(xx, v_length(xx)-1) == 0)
334 v_length(xx)--;
335 hash_insertint(poly1, xx, c);
336 }
337 l_free(tmplist);
338
339 res1 = mult_poly_schubert(poly1, perm, rank);
340 res = monk(maxvar, res1, rank);
341 free_vec_lincomb(res1);
342
343 if (hash_card(poly) != 0)
344 {
345 poly = mult_poly_schubert(poly, perm, rank);
346 lincomb_add_multiple(res, 1, poly, (freekey_t) v_free1, NULL);
347 }
348 hash_free(poly);
349
350 return res;
351 }
352 #endif
353
354
355 void _mult_ps(void **poly, int n, int maxvar, vector *perm, int rank,
356 hashtab *res);
357
mult_poly_schubert(hashtab * poly,vector * perm,int rank)358 hashtab *mult_poly_schubert(hashtab *poly, vector *perm, int rank)
359 {
360 hash_itr itr;
361 void **p;
362 int i, j, n, maxvar, svlen;
363
364 n = hash_card(poly);
365 if (n == 0)
366 return poly;
367
368 p = (void **) amalloc(2 * n * sizeof(void *));
369 i = 0;
370 maxvar = 0;
371 for (hash_first(poly, itr); hash_good(itr); hash_next(itr))
372 {
373 vector *xx = hash_key(itr);
374 j = v_length(xx);
375 while (j > 0 && v_elem(xx, j-1) == 0)
376 j--;
377 v_length(xx) = j;
378 if (maxvar < j)
379 maxvar = j;
380 p[i++] = hash_key(itr);
381 p[i++] = (void *) (long) hash_intvalue(itr);
382 }
383 claim(i == 2 * n);
384 hash_reset(poly);
385
386 svlen = v_length(perm);
387 v_length(perm) = perm_group(perm);
388 _mult_ps(p, n, maxvar, perm, rank, poly);
389 v_length(perm) = svlen;
390
391 for (i = 0; i < n; i++)
392 v_free(p[2*i]);
393 afree(p);
394
395 return poly;
396 }
397
_mult_ps(void ** poly,int n,int maxvar,vector * perm,int rank,hashtab * res)398 void _mult_ps(void **poly, int n, int maxvar, vector *perm, int rank,
399 hashtab *res)
400 {
401 int i, j, c, lnxx, mv0, mv1;
402 hashtab *res1;
403 void *cc;
404
405 if (maxvar == 0)
406 {
407 vector *w = v_new_copy(perm); /* FIXME: OPTIMIZE! */
408 c = (int) (long) poly[1];
409 if (lincomb_add_element(res, c, w, (freekey_t) v_free1) == 0)
410 v_free(w);
411 return;
412 }
413
414 mv0 = 0;
415 mv1 = 0;
416 j = 0;
417 for (i = 0; i < n; i++)
418 {
419 vector *xx = poly[2*i];
420 lnxx = v_length(xx);
421 if (lnxx < maxvar)
422 {
423 if (mv0 < lnxx)
424 mv0 = lnxx;
425 }
426 else
427 {
428 v_elem(xx, maxvar - 1)--;
429 while (lnxx > 0 && v_elem(xx, lnxx-1) == 0)
430 lnxx--;
431 v_length(xx) = lnxx;
432 if (mv1 < lnxx)
433 mv1 = lnxx;
434 poly[2*i] = poly[2*j];
435 poly[2*j] = xx;
436 cc = poly[2*i+1];
437 poly[2*i+1] = poly[2*j+1];
438 poly[2*j+1] = cc;
439 j++;
440 }
441 }
442
443 res1 = hash_new((cmp_t) v_cmp, (hash_t) v_hash);
444 _mult_ps(poly, j, mv1, perm, rank, res1);
445 _monk_add(maxvar, res1, rank, res);
446 free_vec_lincomb(res1);
447
448 if (j < n)
449 _mult_ps(poly + 2*j, n - j, mv0, perm, rank, res);
450 }
451
452
num_inversions(vector * w)453 int num_inversions(vector *w)
454 {
455 int i, j, n, res;
456 n = v_length(w);
457 res = 0;
458 for (i = 0; i < n-1; i++)
459 for (j = i+1; j < n; j++)
460 if (v_elem(w, i) > v_elem(w, j))
461 res++;
462 return res;
463 }
464
perm_group(vector * w)465 int perm_group(vector *w)
466 {
467 int i = v_length(w);
468 while (i > 1 && v_elem(w, i-1) == i)
469 i--;
470 return i;
471 }
472
mult_schubert(vector * w1,vector * w2,int rank)473 hashtab *mult_schubert(vector *w1, vector *w2, int rank)
474 {
475 hashtab *poly;
476 int svlen1, svlen2, w1len, w2len;
477
478 svlen1 = v_length(w1);
479 v_length(w1) = perm_group(w1);
480 w1len = num_inversions(w1);
481
482 svlen2 = v_length(w2);
483 v_length(w2) = perm_group(w2);
484 w2len = num_inversions(w2);
485
486 if (rank == 0)
487 {
488 rank = (((unsigned) -1) >> 1);
489 }
490 else
491 {
492 /* FIXME: one can say exactly when the product is zero. */
493
494 if (2 * (w1len + w2len) > rank * (rank - 1)
495 || v_length(w1) > rank
496 || v_length(w2) > rank)
497 {
498 v_length(w1) = svlen1;
499 v_length(w2) = svlen2;
500 return hash_new((cmp_t) v_cmp, (hash_t) v_hash);
501 }
502 }
503
504 if (w1len <= w2len)
505 {
506 poly = trans(w1, 0, NULL);
507 poly = mult_poly_schubert(poly, w2, rank);
508 }
509 else
510 {
511 poly = trans(w2, 0, NULL);
512 poly = mult_poly_schubert(poly, w1, rank);
513 }
514
515 v_length(w1) = svlen1;
516 v_length(w2) = svlen2;
517
518 return poly;
519 }
520
all_strings(vector * dimvec)521 list *all_strings(vector *dimvec)
522 {
523 int n, ld, i, j, k;
524 vector *str, *cntvec;
525 list *res;
526
527 ld = v_length(dimvec);
528 cntvec = v_new_zero(ld);
529 n = v_elem(dimvec, ld - 1);
530
531 str = v_new(n);
532 j = 0;
533 for (i = 0; i < ld; i++)
534 {
535 while (j < v_elem(dimvec, i))
536 v_elem(str, j++) = i;
537 }
538
539 res = l_new(); /* FIXME: calculate length first */
540 while (1)
541 {
542 l_append(res, v_new_copy(str));
543
544 j = n-1;
545 v_elem(cntvec, v_elem(str, j))++;
546 while (j > 0 && v_elem(str, j-1) >= v_elem(str, j))
547 {
548 v_elem(cntvec, v_elem(str, --j))++;
549 }
550 if (j == 0)
551 break;
552
553 i = v_elem(str, j-1);
554 v_elem(cntvec, i++)++;
555 while(v_elem(cntvec, i) == 0)
556 i++;
557 v_elem(str, j-1) = i;
558 v_elem(cntvec, i)--;
559
560 for (i = 0; i < ld; i++)
561 {
562 for (k = 0; k < v_elem(cntvec, i); k++)
563 v_elem(str, j++) = i;
564 v_elem(cntvec, i) = 0;
565 }
566 }
567
568 v_free(cntvec);
569 v_free(str);
570 return res;
571 }
572
all_perms(int n)573 list *all_perms(int n)
574 {
575 list *res;
576 vector *dimvec = v_new(n);
577 int i;
578 for (i = 0; i < n; i++)
579 v_elem(dimvec, i) = i + 1;
580 res = all_strings(dimvec);
581 v_free(dimvec);
582 return res;
583 }
584
string2perm(vector * str)585 vector *string2perm(vector *str)
586 {
587 int N, n, i, j;
588 vector *dimvec, *perm;
589
590 n = v_length(str);
591
592 N = 0;
593 for (i = 0; i < n; i++)
594 if (N < v_elem(str, i))
595 N = v_elem(str, i);
596 N++;
597
598 dimvec = v_new_zero(N);
599 for (i = 0; i < n; i++)
600 v_elem(dimvec, v_elem(str, i))++;
601 for (i = 1; i < N; i++)
602 v_elem(dimvec, i) += v_elem(dimvec, i-1);
603
604 perm = v_new(n);
605
606 for (i = n-1; i >= 0; i--)
607 {
608 j = v_elem(str, i);
609 v_elem(perm, --v_elem(dimvec, j)) = i;
610 }
611
612 v_free(dimvec);
613 return perm;
614 }
615
str2dimvec(vector * str)616 vector *str2dimvec(vector *str)
617 {
618 int i, n;
619 vector *res;
620 n = 0;
621 for (i = 0; i < v_length(str); i++)
622 if (n <= v_elem(str, i))
623 n = v_elem(str, i) + 1;
624 res = v_new_zero(n);
625 for (i = 0; i < v_length(str); i++)
626 v_elem(res, v_elem(str, i))++;
627 for (i = 1; i < n; i++)
628 v_elem(res, i) += v_elem(res, i-1);
629 return res;
630 }
631
perm2string(vector * perm,vector * dimvec)632 vector *perm2string(vector *perm, vector *dimvec)
633 {
634 int n, i, j;
635 vector *res;
636
637 n = v_elem(dimvec, v_length(dimvec) - 1);
638 res = v_new(n);
639 j = 0;
640 for (i = 0; i < v_length(dimvec); i++)
641 while (j < v_elem(dimvec, i))
642 {
643 int wj = (j < v_length(perm)) ? v_elem(perm, j) : j+1;
644 v_elem(res, wj - 1) = i;
645 j++;
646 }
647
648 return res;
649 }
650
mult_str_schubert(vector * str1,vector * str2)651 hashtab *mult_str_schubert(vector *str1, vector *str2)
652 {
653 vector *dv1, *dv2, *w1, *w2;
654 hashtab *s, *res;
655 hash_itr itr;
656 int i;
657
658 if (v_length(str1) != v_length(str2))
659 return NULL;
660 dv1 = str2dimvec(str1);
661 dv2 = str2dimvec(str2);
662
663 if (v_cmp(dv1, dv2) != 0)
664 {
665 v_free(dv1);
666 v_free(dv2);
667 return NULL;
668 }
669 v_free(dv2);
670
671 w1 = string2perm(str1);
672 w2 = string2perm(str2);
673 /* fixme */
674 for (i = 0; i < v_length(w1); i++)
675 {
676 v_elem(w1, i)++;
677 v_elem(w2, i)++;
678 }
679 s = mult_schubert(w1, w2, v_length(w1));
680 v_free(w1);
681 v_free(w2);
682
683 res = hash_new((cmp_t) v_cmp, (hash_t) v_hash);
684 for (hash_first(s, itr); hash_good(itr); hash_next(itr))
685 {
686 vector *str = perm2string(hash_key(itr), dv1);
687 hash_insertint(res, str, hash_intvalue(itr));
688 }
689
690 free_vec_lincomb(s);
691 v_free(dv1);
692
693 return res;
694 }
695
696